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(ABSTRACT) 


Reliable estimation of thermal properties is extremely important in the utilization 
of new advanced materials, such as composite materials. The accuracy of these estimates 
can be increased if the experiments are designed carefully. The objectives of this study 
are to design optimal experiments to be used in the prediction of these thermal properties 
and to then utilize these designs in the development of an estimation procedure to 
determine the effective thermal properties (thermal conductivity and volumetric heat 
capacity). 

The experiments were optimized by choosing experimental parameters that 
maximize the temperature derivatives with respect to all of the unknown thermal 
properties. This procedure has the effect of minimizing the confidence intervals of the 
resulting thermal property estimates. Both one-dimensional and two-dimensional 
experimental designs were optimized. A heat flux boundary condition is required in both 
analyses for the simultaneous estimation of the thermal properties. For the one- 
dimensional experiment, the parameters optimized were the heating time of the applied 
heat flux, the temperature sensor location, and the experimental time. In addition to these 
parameters, the optimal location of the heat flux was also determined for the two- 



dimensional experiments. 

Utilizing the optimal one-dimensional experiment, the effective thermal 
conductivity perpendicular to the fibers and the effective volumetric heat capacity were 
then estimated for an IM7-Bismaleimide composite material. The estimation procedure 
used is based on the minimization of a least squares function which incorporates both 
calculated and measured temperatures and allows for the parameters to be estimated 
simultaneously. 
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Chapter 1 


Introduction 

A composite material is composed of two or more materials joined together to 
form a new medium with properties superior to those of its individual constituents. There 
are many potential advantages of these materials including higher strength-to-weight 
ratios, better corrosion and wear resistance, and an increased service life over standard 
metals. Because of these improved characteristics, the use of composite materials has 
become quite extensive in the past twenty years, with the most widespread use being in 
the aerospace and aeronautic industries for the design of aircraft structural components. 
For example, composites are used in applications such as aircraft tail sections, wing skins, 
and brake linings. The F-lll horizontal stabilizer was the first flight- worthy composite 
component and in 1986, an all-composite airplane (the Voyager), set a world record in 
nonstop flight around the world, revealing amazing toughness and rigidity against harsh 
environmental conditions. However, the use of composites is not limited to the aerospace 
industry. Composite technology has also gained the attention of the automotive, tooling 
and sporting goods industries. Everything from car bodies and brake linings to tennis 
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rackets, golf clubs, bicycles, and fishing rods have been successfully manufactured from 
composite materials. 

Composites are typically classified according to their reinforcement forms; these 
include particulate, fiber, laminar, flake, and filled/skeletal (Vinson and Sierakowski, 
1987). Fiber-reinforced composites can be further classified as continuous or 
discontinuous. The major types of reinforcing fibers used in composites include glass, 
carbon/graphite, organic, boron, silicon carbide and ceramic fibers, while the major matrix 
resins consist of epoxy, polyimide, polyester, and thermoplastic, with epoxy resins being 
the most versatile of the commercially available matrices. The composite materials 
focused on in this study consist of continuous carbon fiber-epoxy matrix combinations. 

With the increased use of composite materials in aerospace structures and other 
applications, it is important that the properties of these advanced materials be known for 
design purposes. Many studies on the mechanical properties of composites have been 
conducted; however, limited analyses have been made regarding the thermal properties. 
Knowledge of the thermal properties becomes important when the composite is subjected 
to a non-isothermal environment which creates thermal loads on the component. These 
thermal loads induce temperature variations within the structure, which in turn results in 
the development of thermal stresses and possible structural failure. In order to accurately 
predict these thermal stresses and prevent component damage, the temperature response 
of the structure must first be known. However, to determine this response, the thermal 
properties of the composite sample, which can be thermally or directionally dependent, 
are required. The prediction of these thermal properties has provided the motivation for 
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this study. This information will then aid designers in estimating thermal stresses existing 
in a structural component and in turn, allow them to prevent component failure. 

1.1 Goals and Objectives 

The main goal of this research is to predict the thermal properties of composite 
materials. This prediction requires temperature measurements, and therefore, experiments 
must be conducted. The overall objectives of this study are to 

1) develop optimal experimental designs to be used in the prediction of these thermal 
properties 

and 

2) utilize these optimal designs in the development of an estimation procedure to 
determine the effective thermal properties, namely the thermal conductivity and 
volumetric heat capacity. 

Optimal experiments were designed for both isotropic and anisotropic composite materials 
by selecting optimal experimental parameters that maximize the sensitivity of the 
temperature response with respect to changes in the unknown thermal properties. An 
isotropic material has identical properties in every direction while materials exhibiting 
directional characteristics are called anisotropic. For the anisotropic composite material, 
the effective thermal conductivity both parallel and perpendicular to the fiber axis 
direction can be estimated. This optimization procedure was performed because it 
increases the accuracy in the resulting thermal property estimates by minimizing the 
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confidence intervals of the estimated parameters. 

The experimental designs that were optimized not only depend on the boundary 
conditions used, but also on what variability is permitted. An imposed heat flux at one 
boundary, resulting in conductive heat transfer through the composite sample, is required 
in the design to allow for the simultaneous estimation of the thermal properties. 
Therefore, optimal experimental parameters, such as the duration of the applied heat flux, 
should be determined. The optimal experimental parameters determined for the isotropic 
case include the heating time, sensor location, and experimental duration. For the 
anisotropic case, two different experimental designs were used. Both designs had a 
uniform heat flux applied over a portion of one boundary. However, this portion varied 
for the two configurations. Therefore, in addition to the parameters optimized for the 
one-dimensional case, the optimal position of the heat flux was also found in the two- 
dimensional analysis. 

Utilizing the optimal experimental design determined for the isotropic composite 
material, the effective thermal conductivity perpendicular to the fiber axis and the 
effective volumetric heat capacity were then estimated for a composite consisting of 
continuous IM7 graphite fibers and a Bismaleimide (5260) epoxy matrix. Note that this 
is actually an anisotropic composite material; however, since the thermal conductivity is 
only estimated in one direction, this is equivalent to using an isotropic material. The 
estimation procedure used in this investigation was the Gauss linearization method and 
is based on the minimization of a least-squares function, containing experimental and 
calculated temperatures, with respect to the unknown thermal properties. This method not 
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only allows for the effective thermal conductivity and effective volumetric heat capacity 
to be estimated simultaneously, but also enables validation of the transient heat 
conduction equation. 
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Chapter 2 


Literature Review 

2.1 Determination of Thermal Properties of Composite Materials 

This chapter summarizes the present state of knowledge pertaining to the 
estimation of thermal properties of composite materials. Due to their anisotropic nature, 
the estimation of the thermal properties of composites has proved to be a challenging task. 
This estimation problem is further complicated because a composite consists of at least 
two different materials, each with different thermal properties. Many methods, both 
experimental and analytical, have been proposed for estimating these properties with the 
thermal conductivity being most frequently estimated. In the following two sections, 
these estimation techniques are reviewed, describing the methods and procedures used. 
The experimental techniques utilized include both steady-state and transient heat 
conduction processes, while the analytical methods estimate the effective thermal 
properties using proposed mathematical models. These models assume prior knowledge 
of the thermal properties of the fiber and matrix themselves, along with the void fraction 
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of the fibers. The third section describes a minimization procedure based on the Gauss 
method used to estimate the thermal properties. The advantage of this procedure over 
previous techniques is that it allows thermal properties, such as thermal conductivity and 
volumetric heat capacity, to be estimated simultaneously. The thermal properties are 
found by minimizing an objective function containing calculated and measured 
temperatures. The last section discusses optimal experimental designs to be used with this 
minimization procedure which provide more accurate parameter estimates. Optimal 
experimental parameters to be used in these designs are found by maximizing the 
sensitivity of the temperature response with respect to changes in the thermal properties. 

2.1.1 Experimental Determination of the Thermal Properties of Composite Materials 
Experimental methods have been one of the main areas for determining the 
thermal properties of composite materials. These methods can be classified as either 
steady-state or transient Ziebland (1977) described some steady-state experiments used 
to calculate the thermal conductivity that used both absolute measurements, where the 
thermal conductivity is determined directly from the measured quantities, and relative 
methods, in which the thermal conductivity is determined by reference to a substance of 
known thermal conductivity. The absolute methods are accurate but require expensive 
instrumentation and are generally time consuming and thus, expensive. One steady-state, 
absolute technique frequently used is the guarded hot-plate method. In this method, the 
specimen is heated by a hot metal plate attached to it and the resulting temperature is 
measured at the interface to estimate the thermal conductivity (Ziebland, 1977). Although 
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this method is quite accurate, substantial time is required to reach steady-state; therefore, 
the experiment is both expensive and time consuming. 

Dickson (1973) has also described a simple steady-state method for measuring the 
thermal conductivity of insulation materials using heat flow sensors. This method 
requires the measurement of a heat flux and the temperature difference across a test 
specimen of known thickness. Penn, et al. (1986) extended this method to composite 
materials and developed a thermal conductivity measuring apparatus that uses heat flow 
sensors. This steady-state device used smaller sample sizes and as a result, reached 
thermal equilibrium in only a few hours. In addition, Harris, et al. (1982) used a two 
plate apparatus to experimentally determine the thermal conductivities of Kevlar 49 fibers 
in directions parallel and perpendicular to their lengths as functions of temperature, while 
Havis, et al. (1989) experimentally investigated the effect of fiber direction on the 
effective thermal conductivity of fibrous composite materials. 

The evaluation of thermal conductivity from steady-state experiments is 
mathematically simple but frequently lengthy; it was for this reason that transient methods 
were developed. One transient method used to determine the thermal diffusivity, heat 
capacity, and thermal conductivity of materials is the laser-flash method which was first 
introduced by Parker, et al. (1961). In this method, the front face of a small sample is 
subjected to a short, radiant energy pulse. The resulting temperature rise on the rear 
surface of the sample is measured and the thermal diffusivity is then determined from the 
time required for the back surface to reach one half of the maximum temperature rise. 
This can be mathematically expressed as 
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« = El ( 2 . 1 ) 

t 

where K is the constant corresponding to one-half of the maximum temperature rise, L 
is the sample thickness, and t m is the time taken for the back surface to reach one-half 
of the maximum temperature rise. The heat capacity is found from the maximum 
temperature rise of the specimen, and the thermal conductivity is then calculated from the 
product of the thermal diffusivity, heat capacity, and density (k=apc p ). The advantage 
of this technique over steady-state methods is that smaller sample sizes and shorter 
experimental durations could be used. Taylor, et al. (1985) studied the applicability of 
the laser-flash technique for measuring the thermal diffusivity of fiber-reinforced 
composites and found that the technique is appropriate for examining the transient heat 
flow in these materials. 

Lee and Taylor (1975) used the laser-flash method along with an absolute method 
to directly measure the thermal diffusivity of graphite/carbon fiber in unidirectionally 
fiber-reinforced composites. The thermal diffusivity of graphite fiber-reinforced 
composites (Morganite II and Thronal 50 S) was also calculated from the effective 
thermal conductivity of composite samples measured by an absolute method. Taylor and 
Kelsic (1986) also used the laser-flash method to measure the thermal diffusivity of 
unidirectional fiber-reinforced composites. They then investigated the effects of the 
thermal conductivity ratio, fiber fraction, fiber orientation, and specimen length on the 
thermal diffusivity. Their results indicated that the fiber-matrix thermal conductivity ratio 
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was the major factor governing the thermal behavior followed by the fiber volume 
fraction. In addition, the thermal diffusivity of both silica and carbon fiber-phenolic resin 
composites was measured as a function of temperature using the laser-flash technique 
(Mottram and Taylor, 1987a). This work was extended (1987b) and the effective thermal 
conductivity parallel and perpendicular to the fiber axis was calculated using specific heat 
and density data. 

The composite method was used by Brennan, et al. (1982) to measure the thermal 
conductivity and diffusivity of silicon carbide fibers. This method consists of measuring 
the thermal diffusivity of the composite and the matrix itself (without the fibers) using 
the laser-flash technique. From the definition of thermal diffusivity and the Rule-of- 
Mixtures (discussed in the next section), the thermal properties of the fiber can then be 
determined. It was found that the accuracy of the thermal conductivity values determined 
for the fibers could be increased by using a matrix material with a thermal conductivity 
as close as possible to that of the fibers. Furthermore, for this method to yield reliable 
data, it is essential that the scale of the microstructure and the size of the composite 
sample behave as a continuum in its transient response (Brennan, et al., 1982). 

The laser-flash method also served as the basis for the techniques developed by 
Welsh, et al. (1987, 1990). In these studies, a pulsed heat flux was imposed on the 
surface of a material and the resulting thermal response at the same surface was then 
recorded. This method differs from the traditional laser-flash method in that the 
temperature response is observed at the heated surface rather than on the surface opposite 
to the flux. One disadvantage of this method is that the heat capacity and thermal 
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conductivity cannot be estimated independently, only the thermal diffusivity can be 
determined. 

In addition, Fukai, et al. (1991) also conducted transient experiments using a 
periodic hot-wire heating method to simultaneously estimate the thermal conductivity and 
diffusivity. In this method, the thermal conductivity and diffusivity were determined from 
the amplitude and phase lag of the temperature response. The calculated properties agree 
well with those measured by conventional methods. Beck and Al-Araji (1974) also used 
a transient experiment to estimate thermal conductivity and volumetric heat capacity 
independently. 

2.1.2 Mathematical Determination of the Thermal Properties of Composite Materials 

Mathematical models that are functions of the components of a composite have 
also been used to determine the effective thermal properties, particularly thermal 
conductivity. These models are based on the original theories by Maxwell and Rayleigh 
(Hasselman and Johnson, 1987), with the effective properties being direct functions of the 
thermal properties of the constituents, namely the fiber and the matrix. Therefore, it is 
assumed that the thermal properties of the matrix and fiber are known, along with the 
void fraction of the fibers. Typically, results of the mathematical model approach are 
expressed as the ratio of the effective conductivity of the composite to the matrix 
conductivity. This ratio depends on the ratio of the volume of the fiber to the total 
volume and the fiber-matrix conductivity ratio (Han and Cosner, 1981). Hasselman 
(1987) also found that if an interfacial thermal barrier resistance was present in a 
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composite system, the effective thermal conductivity not only depends on the volume 
fraction of the fibers but also on the fiber size. 

One mathematical model, known as the Rule-of-Mixtures, to describe the effective 
thermal conductivity (k eff ) of a composite with heat flow parallel to the axis of the fiber 
is given by 

Kff ~ + (1 ( 2 . 2 ) 


where k f is the thermal conductivity of the fibers, k m is the thermal conductivity of the 
matrix, and V f is the fiber volume fraction. 

A unit-cell approach was presented by Ziebland (1977) to describe the thermal 
conductivity of a composite perpendicular to the fiber axis; this can be mathematically 
expressed as 


k ~ = 

eff 


k k f 

m f 


KVj + a - v f )k f 


(2.3) 


The Rule-of-Mixtures has also been used to calculate the effective thermal 
diffusivity (Taylor and Kelsic, 1986). 


a 


eff 


VA + v'A 

y/pc), + v„(pc) n 


(2.4) 


Here, V m is the matrix volume fraction and (pc)y and (pc) m are the volumetric heat 
capacities of the fiber and matrix, respectively. 

As indicated by Progelhof, et al. (1976), none of the correlations developed 
accurately predict the thermal properties of all types of composites. A review of 
additional models used to predict the thermal conductivity of composite systems is given 
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by Progelhof, et al. (1976). Beran and Silnutzer (1971) presented upper and lower bounds 
for the effective thermal conductivity of a fiber-reinforced composite in terms of volume 
fractions and a geometric factor. They found that the effective thermal conductivity could 
be significantly increased by changing the packing geometry. 

In addition to the analytical models used to estimate thermal properties, numerical 
methods have also been incorporated. Havis, et al. (1989) developed a numerical model 
using the finite difference method that calculated the effective thermal conductivity of 
aligned fiber composites when the fiber to matrix conductivity ratio was greater than one. 
James and Harrison (1992) extended this finite difference method to enable the calculation 
of the temperature distribution and effective thermal conductivity in composite materials 
made from anisotropic materials. The standard finite difference equations were modified 
on a node-by-node basis to take into account anisotropy by local re-orientation of the grid. 
A finite difference method has also been used by James and Keen (1985) to calculate the 
thermal conductivity of uniaxial fiber composites. The effective thermal conductivity was 
then found from the fiber-matrix ratio for a range of fiber volume fractions. This finite 
difference approach was modified by James, et al. (1987) to calculate the transverse 
thermal conductivity of continuous fiber composites in which the fibers can be at any 
angle to the faces of the sample. 

In addition to the finite difference approach, finite elements has also been used to 
predict thermal properties. A finite element analysis of a unit-cell approach was used by 
Han and Cosner (1981) to measure the effective thermal conductivity of fibrous 
composites for two different geometrical arrangements of the fibers, rectangular and 
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staggered. Their analysis assumed prior knowledge of the geometry and thermal 
conductivities of the composite constituents. Veyret, et al. (1993) also used a finite 
element formulation to determine the effective thermal conductivity of a composite 
material using the Laplace equation. 

Other methods have also been used to determine thermal properties. One such 
method is based on the analogy between the response of a unidirectional composite to 
longitudinal shear loading and to transverse heat transfer (Springer and Tsai, 1967). In 
this approach, the thermal conductivities of unidirectional composites were predicted by 
replacing the composite stiffness with the thermal conductivity and the shear modulus 
ratio with the thermal conductivity ratio of the components in the numerical solutions 
obtained for the shear loading problem. Ishikawa (1980) used a method that was 
equivalent to that used by Springer and Tsai. His method was again based on the 
longitudinal shear problem and measured the thermal conductivities of unidirectional, 
carbon-epoxy composite systems using an apparatus based on the infra-red radiation 
method. These analytical results were obtained using a Fourier series analysis and 
required knowledge of the thermal conductivity of the matrix and the fiber volume 
fraction. 

Another technique presented by Behrens (1968) used the method of long waves 
to obtain the average thermal conductivity. By calculating the thermal waves damping 
coefficients in the principal directions of the medium, Behrens was able to develop 
explicit expressions for the average thermal conductivity. In addition, Mottram (1992) 
developed design charts to estimate the effective longitudinal and transverse thermal 
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conductivities of continuous composites using only the fiber and matrix properties. 


2.2 Minimization Methods Used for the Estimation of Thermal Properties 

An alternate procedure for estimating the thermal properties of composite materials 
is to use a minimization technique. One minimization technique frequently used is the 
Gauss linearization method. This is an iterative procedure that involves the minimization 
of the least squares function. Beck (1963) was the first to use this minimization 
procedure to estimate thermal properties, namely thermal diffusivity. 

2.2.1 Gauss Linearization Method 

The Gauss Linearization method, which is based on the least squares function, is 
one of the more popular estimation methods used. This method not only allows for the 
simultaneous estimation of the thermal properties, but also enables validation of the 
transient heat conduction equation. A least squares function, as given by Beck and 
Arnold (1977), is 

S - [Y - T(g)] T [Y - 21®] (2.5) 

where Y is a vector containing measured temperatures, T(Q) is a vector containing 
calculated temperatures, and £ is the true parameter vector. Here, the thermal properties 
are found by minimizing the square of the difference between the measured temperatures 
and the calculated temperatures. For continuous, transient temperature measurements, the 
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sum of squares function is minimized with respect to the parameters using the Taylor 
series approach. This is done by differentiating S with respect to setting the resulting 
equation equal to zero, and then solving for b, the estimated parameter vector for ]3. This 
method, as described by Beck and Arnold (1977), is one of the simplest and most 
effective methods for seeking minima which are reasonably well-defined provided that 
the initial estimates are in the general region of the minimum. However, as explained by 
Box and Kanemasu (1972), if poor initial estimates for the parameters are used or severe 
non-linearity in the model exists, this method may cause large oscillations to occur from 
one iteration to another which leads to non-convergence of the estimates. In an attempt 
to improve the Gauss method, Box and Kanemasu (1972) modified it by changing the 
step size used in seeking the minimum. However, this method still did not include a 
check that the sum of squares function, S , decreased from iteration to iteration. Bard 
(1970) modified the Box-Kanemasu method to include this check; if the function was not 
decreasing, the step size is reduced by one-half. 

The Gauss estimation procedure was used by Beck when he determined the 
thermal conductivity and specific heat of nickel simultaneously from transient temperature 
measurements (1966a) and the thermal contact conductance for both steady-state and 
transient conditions with a periodic contact (1988). Scott and Beck (1992a) also used this 
method to simultaneously estimate the thermal conductivity and volumetric heat capacity 
of carbon composites as functions of temperature and fiber orientation. They found that 
the thermal properties increased with temperature over the range studied and different 
stacking orientations resulted in significantly different thermal conductivity values. This 
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method was also used by Scott and Beck (1992b) to develop an estimation methodology 
for thermoset composite materials during curing, and by Xu and Bao (1990) to 
simultaneously estimate thermal conductivity and diffusivity. 

Loh and Beck (1991) performed a two-dimensional analysis using this estimation 
procedure to simultaneously determine the effective thermal conductivities of anisotropic 
thermoset carbon composites parallel and perpendicular to the fiber axis. They found that 
the conductivity parallel to the fibers is about seven times higher than transverse to the 
fibers. In addition, Jurkowski, et al. (1992) used this method to simultaneously estimate 
the thermal conductivity and thermal contact resistance, as did Gamier, et al. (1992) to 
simultaneously estimate thermal conductivity and volumetric heat capacity without 
internal temperature measurements. Instead, temperature measurements were made with 
thin resistance thermometers and thermocouples. Using finite differences to describe the 
heat transfer model, Pfahl and Mitchel (1970) used this minimization technique to 
estimate six thermal properties of a charring carbon-phenolic material. The calculated 
property values were shown to be in good agreement with values from conventional tests. 

2.3 Optimal Experimental Designs 

Reliable estimation of thermal properties is extremely important in the utilization 
of composite materials. The accuracy of these estimates can be increased if the 
experiments are designed carefully. A carefully designed experiment is one in which 
there is minimum correlation between the estimated properties, as well as maximum 
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sensitivity of the measured experimental variables to changes in the properties being 
estimated (Beck and Arnold, 1977). To create such optimal experimental designs, optimal 
experimental parameters should first be determined. Many criteria have been proposed 
for the design of optimal experiments. Beck and Arnold (1977) have listed some of these 
criteria that are all in terms of the product of the sensitivity coefficients and their 
transpose (X T X). These coefficients are the derivative of temperature with respect to the 
parameters being estimated. The proposed criteria are (1) maximization of the 
determinant of X T X, (2) maximization of the minimum eigenvalue of X T X, and (3) 
maximization of the trace of X T X. The first method was chosen in this study because it 
has the effect of minimizing the confidence intervals of the resulting estimates. 

This optimization method was used by Beck to determine the optimal experiments 
for the simultaneous estimation of thermal conductivity and specific heat (1969) and to 
determine the optimal transient experimental design for estimating the thermal contact 
conductance (1966b). Taktak, et al. (1991) also used this technique to determine the 
optimal heating time of an applied heat flux, optimal number of temperature sensors, and 
the optimal temperature sensor location for the estimation of thermal conductivity and 
volumetric heat capacity of a semi-infinite and a finite thickness composite material. 

As explained, several methods for estimating the thermal properties of composite 
materials have been proposed. These include both experimental methods and the use of 
mathematical models. The procedure used in this study to estimate the thermal properties 
is a modification of the Gauss Linearization method discussed in Section 2.2.1. This 
method was chosen because it allows for the effective thermal conductivity and effective 
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volumetric heat capacity to be estimated simultaneously. Also, when using this technique, 
optimal experiments can be designed resulting in more accurate parameter estimates. 
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Chapter 3 


Theoretical Analysis 

In this chapter, the theoretical development used to determine the optimal 
experimental designs for both isotropic and anisotropic composite materials is presented. 
The minimization procedure used to estimate the effective thermal conductivity 
perpendicular to the fiber axis and the effective volumetric heat capacity of a carbon 
fiber-epoxy matrix composite is also discussed. Recall that this estimation procedure 
requires both experimental and calculated temperatures. In this study, both exact 
analytical temperature solutions and numerical temperature solutions were obtained, with 
the two results being compared to determine the accuracy of the numerical results. The 
numerical solutions were calculated using a finite element program called Engineering 
Analysis Language (EAL, Whetstone, 1983). This finite element software was utilized 
because of the need for future analyses of complex structures, typical in aerospace 
components, for which exact solutions are either complicated or unavailable. 

The first section of this chapter focuses on the mathematical models used to 
describe one-dimensional (isotropic) and two-dimensional (anisotropic) heat conduction 
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processes. The second section describes the mathematical details of the parameter 
estimation technique used in both the exact and numerical analyses to estimate the 
thermal properties. Note that in both cases, the thermal properties estimated were the 
effective properties of the composite, not of the individual fiber and matrix components. 
In the final section, the mathematical criterion used to design optimal experiments, 
resulting in greater accuracy of the thermal properties, is discussed. 

3.1 Mathematical Models Used in Estimating the Thermal Properties of 

Composite Materials 

The formulation of a mathematical model is based on the experimental system 
being analyzed. In this investigation, formulating mathematical models, either exact or 
numerical, to describe the conductive heat transfer occurring within the composite sample 
will allow for the temperature distribution to be calculated. This distribution is required 
for the estimation of the thermal properties. As mentioned, both one-dimensional and 
two-dimensional heat conduction analyses have been conducted. The mathematical 
formulation behind both are defined in the following two subsections. 

3.1.1 One-Dimensional Analysis - Isotropic Composite Material 

For the isotropic situation, one-dimensional heat conduction through a carbon- 
epoxy composite was investigated. Note that this isotropic situation is equivalent to 
analyzing the properties in one direction of an anisotropic composite, as was the case in 
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this study. The samples used consisted of a thin, flat disk with an aspect ratio such that 
the two-dimensional effects at the edges can be ignored. One plane boundary had an 
imposed heat flux perpendicular to the fiber axis, and a known, constant temperature 
existed at the second boundary, as shown in Fig. 3.1. Since composite materials tend to 
have low thermal conductivities in directions perpendicular to the fibers, this isothermal 
boundary condition is readily available. The heat flux boundary condition is required for 
the independent estimation of the thermal properties. This requirement occurs because 
this type of boundary condition introduces a new equation into the model which contains 
only the thermal conductivity and not the volumetric heat capacity. This equation is 
known as Fourier’s Law and is given by 



where q x is the applied heat flux. If this boundary condition was not used, and instead, 
a constant temperature or insulated condition was used, then only the thermal diffusivity 
(k/pc p ) could be estimated. 

The formulation to describe this problem can be found from an energy balance and 
is expressed as 


3 (lr dT 1 

3 * 1 ^^ 



0 < x < L x 


t > 0 


(3.2) 


where T is temperature, k x ^ ff and x are the effective thermal conductivity and position, 
respectively, in the direction of heat transfer, C eff is the effective volumetric heat capacity 


22 




Figure 3.1. Experimental Set-Up for the Estimation of the Effective Thermal 

Conductivity and Volumetric Heat Capacity of an Isotropic Material. 
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(or product of density and specific heat), and t is time. The heat flux and constant 
temperature boundary and initial conditions can be described as 
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where t h is the time that the heat flux is applied to the sample. After this time, the 
boundary condition becomes insulated, as seen by Eqs. (3.3a, b). The heat flux, q x , the 
temperature at x = L x (T ox ), and the initial temperature, 7J, are assumed to be known 
without errors. Note that two solutions were required for this analytical problem; one 
while the heat flux was applied and one after the duration of the heat flux. Also, since 
the experiments were conducted at room temperature, it was assumed that the temperature 
at x = L x was equal to the initial temperature; i.e. T ox - 7j. Using these assumptions, the 
exact solutions to describe the temperature distributions were obtained using Green’s 
function (Beck, et al., 1992). The Green’s function required for this solution is given by 
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where p n is an eigenvalue represented by 
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(Beck, et al., 1992). Using Eq. (3.6), the one-dimensional temperature distribution was 
solved for, resulting in the following: 
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for t > L. 


The temperature solution was also obtained numerically from the finite element 
software, EAL, using an implicit transient analysis. In EAL, the weighted residual 
method is used to derive the implicit time integration equations. During each time step, 
the temperature vector is approximated by 


(C + A tK)T i+1 = (C - A tK)T. + FAt + FAt 2 (3.10) 

where T ( is the temperature vector at time T i+1 is the temperature vector at time t i+J , At 
is the time step size, C is the capacitance matrix, K is the stiffness matrix, and F is the 
matrix containing the boundaiy conditions. (Whetstone, 1983). Again, a numerical 
approach was utilized for the future need to analyze complex structures which do not 
have exact solutions available. 


3.1.2 Two-Dimensional Analysis - Anisotropic Composite Material 

The two-dimensional analysis is similar to the one-dimensional analysis, only now. 
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two-dimensional heat conduction through an anisotropic composite sample is considered. 
Two different experimental configurations were used in this analysis. The first consists 
of an imposed heat flux perpendicular to the fiber axis over a portion of one boundary 
(with the remainder of the boundary insulated) and known constant temperatures at the 
remaining three boundaries, as shown in Fig. 3.2. The second configuration also has a 
heat flux imposed over a portion of one boundary, only now, the boundary opposite to 
the heat flux is maintained at a constant temperature, while the remaining two boundaries 
are insulated, as shown in Fig. 3.3. For both experimental assemblies, the heat flux 
boundary condition will allow for the determination of thermal conductivity in two 
directions. However, the actual estimation of these thermal conductivities will not be 
performed in this study; only the experimental designs required for this estimation process 
will be analyzed (Section 3.3). 

The temperature distribution within the material for both configurations can be 
determined from conservation of energy 


a (, dT\ d(, dT) „ dT 

dx{ x - ,r te J dy [ y - er ty J eir di 


0 < x < L. 0 < y < L t > 0 (3.11) 

y 


where, in this case, k y , eff and y are the effective thermal conductivity and position, 
respectively, perpendicular to the direction of heat transfer. The temperature solutions 
obtained for both configurations are discussed in the following two subsections. 
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Figure 3.2. Experimental Set-up Used for Configuration 1 in the Estimation of the 
Effective Thermal Conductivities in Two Orthogonal Planes. 



Figure 3.3. Experimental Set-up Used for Configuration 2 in the Estimation of the 
Effective Thermal Conductivities in Two Orthogonal Planes. 
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3. 1.2.1 Configuration 1 - Isothermal Boundary Conditions 

The heat flux and isothermal boundary conditions and the initial temperature 
condition for Configuration 1 (Fig. 3.2) can be described as 
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where q x is the applied heat flux, T ox , T oyl , and T oy2 are the known temperature boundary 
conditions, L x is the thickness of the plate in the x direction, L y is the thickness of the 
plate in the y direction, is the portion of the plate where the heat flux is imposed, and 
T { is the initial temperature. The specific value for Z^ ; will be found using the 
optimization procedure discussed in Section 3 of this chapter. Note that once again, two 
solutions are required for this analytical case; one while the heat flux is applied and one 
after the duration of the heat flux. Also, since the experiments were again conducted at 
room temperature, it is assumed that T ox = T oyl = T oy2 = T { . Using these assumptions, the 
solutions to describe the temperature distribution within the composite sample were 
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obtained using Green’s functions. For the two-dimensional case, two Green’s functions 
are required for the temperature solution; one for both the x and y direction boundary 
conditions. The Green’s function for the heat transfer along the x axis is provided in Eq. 
(3.6), and the Green’s function along the y axis is given by (Beck, et al., 1992). 
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where is the effective thermal conductivity ratio, (k y _ eJ /k x _ eff ). Using these Green’s 
functions, the temperature solutions for Configuration 1 are represented by 


- T + sin 

f > 
mny 

cos 

V 

1 - cos 

( J \1 

nrnL pl 

OJC r Z-f JLi 

x-eff 11 »-* 


w 

L 

l * p 


V mB J 


[l - exp(-At)] 


for 0 < t < t hy and 


AqL ~ ~ 
T(x,y,t) = T of + ZL 


m - 1 n = 1 

f l ' 


f \ 


(n } 


- 

( Vi 

mny 

cos 

fin* 


1 - cos 

mnL p,l 

L 

L 


L 

V y J 


V x J 


- 

l y P 


ytnB j 


[exp[-A(^-^)j - Qxp(-At)] 


for t > t h , where 


(3.19) 


(3.20) 


A = 


m 


V k 


B 2 k 

y-eff i Fn \- e ff 


l 2 c 


LlC 


<# ) 


(3.21) 


B = mViiic + p 


xyxy t'n 


(3.22) 


29 



and L xy is the ratio of the composite dimensions (L/L y ). 


3. 1.2.2 Configuration 2 - Isothermal and Insulated Boundary Conditions 

The heat flux and isothermal boundary conditions and the initial temperature 
condition for Configuration 2 (Fig. 3.3) can be described as 
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where L p2 is the portion of the plate where the heat flux is imposed. Again, the specific 
value for L p2 will be found using the optimization procedure discussed in Section 3 of this 
chapter. Due to the different boundary conditions used in the two configurations, L p 2 will 
be different than Lp j (the heat flux position calculated for Configuration 1). Since the 
experiments were again conducted at room temperature, the same assumption was used 
as for Configuration 1; i.e., T ox - T t . The solutions to describe the temperature 
distribution within the composite sample were then obtained using Green’s functions. 
Two Green’s functions are again required for this configuration, one for both the x and 
y direction boundary conditions. The Green’s function along the x axis is provided in Eq. 
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(3.6), and the Green’s function along the y axis is given by (Beck, et aL, 1992) 
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Using these Green’s functions, the temperature solutions for Configuration 2 are 
represented by 
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for t > t h , where A and B are given by Eqs. (3.21) and (3.22), respectively, and C is 
represented by 


C = 


&K-. 


x-eff 


LlC 


‘ff 


(3.32) 


In determining these temperature distributions as functions of time, one should note that 
there are steady state terms which need to be calculated only once since they are time 
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invariant. This is important since these series are slow to converge and require hundreds 
of terms, whereas the time varying terms of the summation converge rather quickly. An 
alternate solution method involves the use of time partitioning (Beck, et al., 1992). In 
this method, the solution is partitioned into two regions and both large-time and small- 
time Green’s functions are used to find the temperature. For example, at early times, the 
solution is the same as that for a semi-infinite body, and therefore, the overall solution 
can be divided up into early and steady state solutions. 

3.2 Minimization Procedure Used in Estimating the Thermal Properties 

The method used to estimate the thermal properties is based on the minimization 
of an objective function with respect to the unknown parameters, effective thermal 
conductivity and effective volumetric heat capacity. This procedure is called the Gauss 
method and allows for the simultaneous estimation of the thermal properties. A 
modification of the Gauss method that allows for nonlinearities in the model to exist is 
the Box-Kanemasu method, which is utilized in this investigation. In this method, the 
objective function used is the least-squares function, S, and is given by 

5 - [Y - 7<g)f [Y - 7(g)] (3.33) 

(Beck and Arnold, 1977). Here, Y is the measured temperature vector, 7(g) is the 
calculated temperature vector found using a transient mathematical model (as given in 
Eqs. (3.8-9), (3.19-20), and (3.30-31)) and the parameter estimates, and g is the exact 
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parameter vector that contains the unknown thermal properties. The objective function, 
S, is minimized with respect to the unknown parameters, ]3. This is done by 
differentiating S with respect to and setting the resulting equation equal to zero, giving 

V = 2[-X T (&] [Y - = 0 (3.34) 

(Beck and Arnold, 1977). Here, the sensitivity coefficient matrix, X(Q) is defined as 

X(S) = [V/W (3.35) 

These coefficients are the derivatives of temperature with respect to the parameters being 
estimated and represent the sensitivity of the temperature response with respect to changes 
in the unknown parameters. In order for the parameters to be estimated simultaneously, 
the determinant of the sensitivity coefficients and their transpose, I X T X | , cannot equal 
zero. That is, any one column of X cannot be expressed as a linear combination of any 
other column. 

Because the heat conduction process in this study is a non-linear problem, the 
estimator, J5, cannot easily be solved for. Therefore, two approximations are used in Eq. 
(3.34) to prevent this difficulty; (1) Replace X(Q) by X(b), where b is an estimate of Jl, 
and (2) Use the first two terms of a Taylor series for T(Q) about b to approximate T(Q) 
(Beck and Arnold, 1977). Using these approximations and implementing an iterative 
scheme, Eq. (3.34) can be solved for b, the estimated parameter vector, resulting in the 
following expression for b (k+1) : 

&(* + D = b <*> + (Y - I™ )] (3.36) 

where 
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j*k) = 


( 3 . 37 ) 


This is known as the Gauss linearization equation Here, k is the iteration number, b (k+1) 
is the new parameter estimate, b (k) is the estimate at the previous iteration, and T(b ) (k> 
contains temperatures calculated using b (k) . 

For a nonlinear problem, Eq. (3.36) is altered and becomes 

£(* +1 > = £(*> + h )a ft<*) (3.38) 

where 

A b (k) = ( Y - r k) )] (3.39) 

and h (k+1) is a scalar interpolation function. To use this nonlinear estimation procedure, 
an initial estimate, b (0) , is required. This estimate is then used to calculate T* oi and 
which are used to obtain the improved parameter vector, b (I) . This procedure continues 
until all parameters in b do not change significantly (Beck and Arnold, 1977). 

Equation (3.38) represents the Box-Kanemasu method which is a modification of 
the Gauss method. In the Box-Kanemasu method, the sum of squares, S, is approximated 
at each iteration by a quadratic function in h. The minimum S is located where the 
derivative of S with respect to h is equal to zero, or at an h value of (Beck and Arnold, 
1977) 

h (kH) = GV[Sf - + IG^aY' (3.40) 

where 

G (k) = [A^Wf (3.41) 
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The parameter for a is initially set equal to one and S a (k) and Sj k) are the values of S at 
a and zero, respectively. If S a (k) is not less than Sj k) , a is reduced by one-half and the 
inequality is checked again. This is a modification over the original Box-Kanemasu 
method. A flow chart illustrating the modified Box-Kanemasu estimation procedure, as 
presented by Beck and Arnold (1977), is shown in Fig. 3.4. Note that if the investigation 
requires a to become less than 0.01, the calculations are terminated. One reason why this 
may occur is that correlation (or linear dependence) between the sensitivity coefficients 
exists, causing the sum of squares function not to have a unique minimum. It is therefore 
very important to calculate and analyze the sensitivity coefficients for possible correlation 
to ensure reliable parameter estimates. 

A parameter estimation program was written using the modified Box-Kanemasu 
method and is called MODBOX ; this program is based on the original program NL1NA, 
by Beck (1993). This program uses sequential in-time estimation to calculate the 
parameters at each time step. The exact mathematical models given in Eqs. (3.8) and 
(3.9) were used in this program as well as the derived sensitivity coefficients, allowing 
for the estimation of the effective thermal conductivity perpendicular to the fibers and the 
effective volumetric heat capacity of a composite consisting of IM7 graphite fibers and 
a Bismaleimide epoxy matrix. The modified Box-Kanemasu method was also 
implemented into EAL where the temperature solution was obtained numerically. Again, 
the same effective thermal properties were estimated. The advantage of this sequential 
estimation technique is that it allows the user to observe the effects of additional data on 
the sequential estimates and study the validity of the proposed mathematical model and 
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Figure 3.4 Flow Chart for the Modified Box-Kanemasu Estimation Procedure 
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experimental design. Ideally, at the conclusion of an experiment, any additional data 
should not affect the parameter estimates. 

3.3 Optimal Experimental Designs Used in Estimating Thermal Properties of 

Composite Materials 

Since the Gauss method requires experimental temperatures of the composite 
system to be measured, the accuracy of the thermal properties estimated can be greatly 
increased if these experiments are designed carefully. To create such optimal 
experimental designs, optimal experimental parameters must first be determined. The 
focus of this section is on the criterion used in obtaining these optimal parameters. For 
the one-dimensional analysis, the experimental design consisted of a thin plate with an 
imposed heat flux applied for a finite duration at one boundary and a known, constant 
temperature at the second boundary. For this design, the optimal experimental parameters 
that were determined are the heating time, temperature sensor location, and total 
experimental time. 

For the two-dimensional heat conduction analysis, two different configurations 
were used, allowing for the effective thermal conductivity in two directions and the 
effective volumetric heat capacity to be determined simultaneously. Both designs had a 
heat flux imposed over a portion of one boundary, with the remainder of the boundary 
insulated. In addition. Configuration 1 had known, constant temperatures at the remaining 
three boundaries, while Configuration 2 had a constant temperature at the boundary 
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opposite to the heat flux, and insulated conditions at the remaining two boundaries. 
Therefore, in addition to the optimal experimental parameters found for the one- 
dimensional case, the optimal position of the heat flux was also determined for both 
configurations used in the two-dimensional case. Note, however, that this optimal heat 
flux location will not be the same for both configurations due to the different boundary 
conditions used. 

3.3.1 Design Criterion Used for Optimal Experimental Designs 

Many criterion have been proposed for the design of optimum experiments. As 
mentioned previously, the sensitivity coefficients indicate the sensitivity of temperature 
to changes in the thermal properties and optimal experiments are those which maximize 
these coefficients for each property. Therefore, the criterion chosen for this analysis is 
the maximization of the determinant (D + ) of X +T X + , which contains the product of the 
dimensionless sensitivity coefficients and their transpose (Beck and Arnold, 1977). This 
criterion is subject to a maximum temperature rise, a fixed number of measurements, and 
the following seven standard statistical assumptions: additive, zero mean, constant 

variance, uncorrelated normal errors with errorless independent variables, and no prior 
information. It is recommended by Beck and Arnold because it has the effect of 
minimizing the confidence intervals of the resulting parameter estimates. Note, it was 
desired to perform the optimization procedure in non-dimensional terms so the results 
could be applicable for any material, not just composite materials. 

For the one-dimensional analysis where two properties are estimated (k x _ eff and Q). 
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I x+ T x + 1 is a 2 x 2 matrix. Therefore, the dimensionless determinant is given as 


D\-d - 


= - (d; 2 y 


du d l2 
dn d } 2 

where d n + , d 12 + , and d 22 + are found from (Beck and Arnold, 1977) 


(3.42a, b) 




i+2 


M t ; 


£ f' ; x,V)x/(n dr 

p* 1 


(3.43) 


In this equation, M is the number of temperature sensors used and t + , t N + , and T max + are 
defined as 


Tl 


(T ma x T 0 J 

QxWK-eff ’ 




(3.44a-c) 


where t N is the total experimental time, T max is the maximum temperature reached between 
the start and end of the experiment, and T ox is the surface (and in this case, initial) 
temperature. It should be noted that this definition of T max + was used to verify previous 
optimal experimental parameter results by Taktak, et al. (1991) and is not the best 
representative choice. 

The integral in Eq. (3.43) was calculated numerically, being approximated by a 
summation. From Eq. (3.43), it is evident that the matrix in Eq. (3.42a) is symmetric; 
i.e., d l2 = d 2} + . This simplifies the problem by decreasing the number of equations that 
must be numerically integrated. 

When extending this analysis to the two-dimensional case, three properties {k x -e& 
k y _ eff and C eff ) can now be estimated simultaneously for both configurations. Therefore, 
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X +r X + I is a 3 x 3 symmetric matrix and the dimensionless determinant is given by 


Kd = 


d\ i d n d\$ 
d\2 d} 2 d^2 

du <^23 ^33 


^ 2 -d = dnid ^ d ^ d > z $) ^12(^12^33 ^13^23) + ^13(^12^23 d ^ dyi ) 


(3.45a,b) 


Again, the values were found from Eq. (3.43), where the integral was calculated 
numerically. To compare both configurations used in the two-dimensional analysis, the 
value for T max + was redefined as the temperature reached at steady-state conditions. This 
is a more accurate choice than the 7’ m<u + selected for the one-dimensional design, used by 
Taktak, et al. (1991), because it represents the true maximum temperature that can be 
attained for the defined problem. 

From Eqs. (3.42) and (3.43), it is seen that the dimensionless sensitivity 
coefficients are required for this optimization procedure; these coefficients are given by 


and 


dT 

c *jr ar 


(3.46) 

(3.47) 

(3.48) 


where X k and Xr are used in the one-dimensional analysis and X^ , X k , 
and Xc g are used in the two-dimensional analysis. 
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3.3.2 One-Dimensional Optimal Experimental Design Formulation 


In performing this optimization procedure, a mathematical model, either exact or 
numerical, is required to represent the experimental process. An exact model for the one- 
dimensional analysis is given by the temperature solutions in Eqs. (3.8) and (3.9). Using 
the following dimensionless variables 
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u = 
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+ _ x-efrh 
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C ej^x C eJ gL? 4x L JK-e ff 

these temperature distributions can be expressed in dimensionless form as 


(3.49a-d) 


T\xj) = 1 - * + - 2£ J_ cos(P n x + ) exp(-pjr) 
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(3.50) 


T\x,t) = -2 £ — cos(P„x*) expc-ftjo - exp[-p*(< *-*,*)]] 

-> P» J 


(3.51) 


for > t h + , where P„ is given by Eq. (3.7). This temperature distribution, which is 
calculated using a dimensionless heating time, t h + , equal to the total experimental time, 
is shown in Fig. 3.5 for several x + locations. 

The dimensionless sensitivity coefficients, X? and , were then found by 


differentiating Eqs. (3.50) and (3.51) with respect to k x . eff and C eff . The effective thermal 
conductivity sensitivity coefficients are given by 
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(3.52) 
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Dimensionless Temperature, 



Dimensionless Time, f 


Figure 3.5. Dimensionless Temperature Distribution (T*) for One-Dimensional Heat 
Conduction Using a Dimensionless Heating Time, t h + , Equal to the Total 
Experimental Time for Several x + (x/L x ) Locations. 



for 0 < t + < t h + , and 


- 4 + o 
p j 


for f > 4 + , while the effective volumetric heat capacity sensitivity coefficients can be 
expressed as 


*c, = -2j^/*exp(-pjr) cos(P n x *) (3.54) 

n-1 

for 0 < f < t h + , and 

*c„ = -2^[(*exp(-(5^*) - (/*- Oexpt-Pfo'-O]] cos(P B x*) (3.55) 

n=l J 

for f > 4 + . 

These dimensionless sensitivity coefficients are then used in the optimization procedure 
to determine the maximum determinant value, as given by Eqs. (3.42) and (3.43), and the 
corresponding optimal experimental parameters. Viewing these coefficients as functions 
of experimental time will also give insight into the experimental design, as will be shown 
in Chapter 5. The sensitivity coefficients in dimensional form are also required in the 
program, MOD BOX, to estimate the effective thermal conductivity and effective 
volumetric heat capacity simultaneously, as shown in Eqs. (3.36) and (3.37). 


3.3.3 Two-Dimensional Optimal Experimental Design Formulation 

The exact model used to describe the temperature distribution for two-dimensional 
heat transfer in an anisotropic composite material is given in Eqs. (3.19) and (3.20) for 
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Configuration 1 and Eqs. (3.30) and (3.31) for Configuration 2. The derived sensitivity 
coefficients for both configurations required for the optimization procedure are discussed 
in the following two subsections. 


3.3.3. 1 Optimal Experimental Design Formulation for Configuration 1 

Using the dimensionless variables already given in Eq. (3.49) along with the 
following non-dimensional variables 


>'-T' i=u 


(3.56a,b) 


y y 

(i corresponds to the configuration number) the temperature distributions for 
Configuration 1 can be expressed in dimensionless form as 
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for 0 < f < 4 + , and 
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• jexp[-if(f + - 4 + )j - exp(-B? + )j (3.58) 

for t + > t h + , where L xy (L/L y ) and k x> (k y . e /k x . eff ) are the dimension and effective thermal 
conductivity ratios, respectively, and B is given in Eq. (3.22). 

The sensitivity coefficients for all three effective parameters, thermal conductivity 
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perpendicular to the fibers (k x . eff ), thermal conductivity parallel to the fibers (k re ^ y and 
volumetric heat capacity (C eff ), were then calculated by differentiating the above 
temperature distribution with respect to each property. The sensitivity coefficients for the 
thermal conductivity perpendicular to the fiber axis are given by 


A °° °° / \ 

sin(m7ty + ) cos(P n Jc + ) [1 - co simrcLp^)] L_ ] 
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for f > t h + where D is equal to 
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(3.61) 


The dimensionless sensitivity coefficients for X,' and X* were also calculated: the 

1 1 y-*ff 1 

solutions for 0 < t + < t h + are given by 
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where £ and D are given in Eqs. (3.22) and (3.61), respectively. 


3.3.3.2 Optimal Experimental Design Formulation for Configuration 2 

Using the dimensionless variables given in Eqs. (3.49) and (3.56), the temperature 
distribution obtained for Configuration 2 (Eqs. (3.30) and (3.31)) can be expressed in 
dimensionless form as 
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for 0 < f < t h + , and 
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for f > r A + . 

Again, the sensitivity coefficients for all three effective parameters, thermal 
conductivity perpendicular to the fibers (k x _ eff ), thermal conductivity parallel to the fibers 
iky.#), and volumetric heat capacity (C eJf ), were then calculated by differentiating the 
above temperature distributions with respect to each property. The sensitivity coefficients 
for the thermal conductivity perpendicular to the fiber axis are given by 
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for t > t h + . The dimensionless sensitivity coefficients for and XC 

^ tf t- 

calculated; the solutions for 0 < f < t h + are given by 
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where B and D are given in Eqs. (3.22) and (3.61), respectively. 

Again, as in the one-dimensional analysis, the dimensionless sensitivity coefficients 
for both configurations will be used in the optimization procedure to determine the 
maximum determinant and the corresponding optimal experimental parameters. 
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Chapter 4 


Experimental Procedures 

This chapter describes the experimental procedure used to estimate the thermal 
properties of a continuous IM7 graphite fiber - Bismaleimide epoxy matrix composite 
material. Although optimal experiments were designed for both one-dimensional 
(isotropic) and two-dimensional (anisotropic) heat conduction, only the one-dimensional 
experiment was conducted, allowing for the effective thermal conductivity perpendicular 
to the fibers and the effective volumetric heat capacity to be estimated simultaneously. 

As discussed previously in Section 3.2, the estimation procedure used in this study 
is the modified Box-Kanemasu method. Recall that when using this method, experimental 
temperatures must be recorded. To estimate the thermal properties independently, the 
experiments must be transient and one of the boundary conditions must be a heat flux 
(Beck and Arnold, 1977). With these required conditions, the experimental assemblies 
were designed accordingly. Discussed next are the experimental set-up and procedure 
utilized. 
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4.1 One-Dimensional Experiment for the Estimation of Thermal Properties 


The experimental assembly for the one-dimensional analysis used to estimate the 
thermal properties of the given composite material consists of a thin composite sample 
subjected to a heat flux perpendicular to the fiber axis at one boundary and a known 
constant temperature at the other boundary. Temperature measurements were then taken 
at the flux boundary and were used in the estimation procedure. These experiments are 
described in detail in the following subsections. 

4.1.1 One-Dimensional Experimental Set-Up 

The experimental design for one-dimensional heat conduction was composed of 
two composite disks of approximately equal size, a resistance heater, eight thermocouples, 
and two copper cylinders. The assembly was symmetrical consisting of, from the center 
to the top, a thin resistance heater, two thermocouples, the composite sample, two 
additional thermocouples, and a copper block. The composite sample was 4.77 cm in 
diameter and 0.678 cm thick. The copper blocks, each with a height of 6.35 cm and a 
diameter of 5.08 cm, were used as heat sinks to attempt to attain the constant temperature 
boundary condition while the resistance heater was used to provide the heat flux boundary 
condition. All of the experiments were conducted at the National Aeronautics and Space 
Administration - Langley Research Center, Aircraft Structures Branch (NASA-LaRC,ASB) 
using the equipment available in their testing lab. All supplies required in the experiment, 
such as the resistance heater and heat sink compound, were also supplied by NASA- 
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LaRC. The carbon fiber-epoxy matrix composites were prepared and the thermocouples 
(Type K) were fabricated by NASA-LaRC personnel. The data acquisition hardware and 
software used in taking temperature, voltage, and current measurements had been 
previously programmed. Therefore, only the assembly of the experimental apparatus 
remained to be completed, with the details given next. 

4. 1 . 1 . 1 Experimental Set-Up Assembly 

The one-dimensional experimental apparatus involved the following procedure: 

1) Measure the thickness and diameter of two composite samples (Samples 1 and 2), 
and the height and diameter of two copper blocks. 

2) Coat one surface of a copper block with a thin layer of silicon heat sink 
compound. Make sure the compound is smooth and evenly distributed. A flat 
edge the width of the sample is useful to apply the coating with. 

3) Place (2) thermocouples on top of the copper block layered with the heat sink. 
The junction should be in the center of the samples and the two wires should be 
parallel to each other and equidistance approximately 0.635 cm from the center of 
the second axis. To keep the thermocouples in place, tape them down either to 
the table or to the sides of the sample. Be sure to number the thermocouples so 
that their position can be recorded (see Fig. 4.1). 

4) Coat one surface of Sample 1 with a thin layer of silicon heat sink compound. 
Again, be sure the compound is smooth and evenly distributed. Carefully place 
the composite sample (coated surface down) on top of the copper block over the 
thermocouples (see Fig. 4.2). Do not slide the composite sample on the block or 
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diameter = 0.0477 m 



T/C 1 


T/C Junction 
Copper Block 


T/C 2 


Figure 4.1. Position of Thermocouples (T/Cs) on Copper Block for the 
One-Dimensional Experimental Design. 



Composite 
Copper Block 


Figure 4.2. Sample 1 Placed on Top of the Copper Block for the 
One-Dimensional Experimental Design. 
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the thermocouples will be moved. Make sure that the contact is good. 

5) Coat the other side of Sample 1 and again, place two thermocouples at this surface 
in the same manner as discussed previously. 

6) Apply silicon grease to the heater on one side and place it on top of Sample 1. 
Be sure the heater is placed symmetrically over the sample so the same magnitude 
of heat flux is being distributed (see Fig. 4.3). 

7) Coat the exposed top surface of the heater with the heat sink compound as 
before and place two more thermocouples at this surface, as described previously. 

8) Coat Sample 2 (of approximately the same thickness as Sample 1) with silicon 
grease and place it on top of the exposed heater surface, over the thermocouples. 

9) Coat the opposite side of Sample 2, place two more thermocouples as before on 
the surface, and finally, place a coated copper block over the composite sample 
(see Fig. 4.4). 

10) Wrap the exposed sides of the composite material with rope insulation. Four 
pieces were used, two on each composite. If any thermocouple wire is exposed, 
tuck it inside of the insulation (see Fig. 4.5). 

11) Carefully place the stacked samples between two plates and apply pressure evenly 
over the surface, taking care not to break the thermocouples. If thermocouples 
break, use less pressure. (Note, pressure was applied through threaded rods which 
ran through the comers of the plates). 
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Figure 4.3. Position of the Heater and Thermocouples (T/Cs) at the Heat Flux 
Boundary Condition for the One-Dimensional Experimental Design. 
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Figure 4.4. Final Assembly of the Experimental Apparatus for the One-Dimensional 
Experimental Design with Eight Thermocouples (T/Cs). 
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Figure 4.5. Insulation Wrapped Around the Heater and Composite Samples 
for the One-Dimensional Experimental Design. 
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4. 1 . 1 . 2 Experimental Procedure 

The experimental procedure consisted of applying a heat flux to a composite 
sample already at steady-state and then measuring the resulting temperatures using the 
data acquisition system. The process used for these experiments is given as follows: 

1) Place the press containing the experimental set-up inside a temperature controlled 
oven. Connect the thermocouple wire to leads leading to a temperature 
compensator and data acquisition system. Heat the oven to the desired ambient 
temperature and allow the instrumented samples to equilibrate. 

2) Activate the data acquisition system. Turn on the heater and simultaneously record 
temperature (mvolts) from thermocouples, and measured voltage and current to the 
heater. Turn the heater off after a pre-determined heating time and continue 
recording measurements until a pre-determined experimental time has elapsed. 

3) If desired, the experiment can be repeated after the samples have again come to 
equilibrium with the oven temperature or the oven temperature can be changed and 
steps 1 and 2 can be repeated. 

In these experiments, temperatures were recorded using eight Type K 
thermocouples at 0.5 second intervals up to a predetermined experimental time. All 
experiments were conducted at room temperature with the heater being applied for a 
predetermined heating time. (These times were determined using the optimal design 
criterion discussed previously). Experiments were conducted using three different voltage 
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inputs to the heater; 4.9V, 6.1V, and 7.3V. These resulted in maximum temperature rises 
of approximately 2°C, 3°C, and 4.5°C, respectively, over the initial temperature. 

As mentioned, a one-dimensional heat conduction process through the composite 
sample was assumed. This assumption may introduce experimental error into the 
problem. However, to verify the validity of this assumption, the resistance to heat 
transfer was calculated both parallel and perpendicular to the direction of heat transfer. 
For one-dimensional heat conduction to be assumed, the parallel resistance should be 
much smaller than the perpendicular resistance, indicating that most of the heat will be 
conducted in one direction. To calculate the resistance in the direction of heat transfer, 
where only conduction through the sample is considered, the following equation was used: 


R 


cond 



(4.1) 


where A x is the cross-sectional area normal to the direction of heat transfer, k x _ eff is the 
effective thermal conductivity parallel to the direction of heat transfer and L x is the 
thickness of the sample. For the resistance perpendicular to the direction of heat transfer, 
both conduction through the insulation material and convection with the air must be 
considered. The convective resistance is given as 


R 


conv 


1 

hA 


(4.2) 


where h is the convective heat transfer coefficient. As an extreme case, this resistance 
was calculated using a h of 10 W/m 2o C. This indicates that heat is lost through the 
insulation by convection to the surrounding air. Since the experiments were conducted 
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at room temperature, this h value is a good approximation for natural convection 
situations, as may be the case in this experiment. A thermal conductivity of 0.05 W/m°C 
was assumed for the insulation material used. 

When performing these calculations, is was found that the resistance parallel to the 
direction of heat transfer is 7.3 °C/W while the resistance normal to this direction is 220 
°C/W. Since the perpendicular resistance is much larger than the parallel resistance, the 
one-dimensional heat conduction assumption is valid. 
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Chapter 5 


Results and Discussion 

This chapter focuses on the results obtained for the optimal experimental design 
procedure for both the one-dimensional and two-dimensional analyses. In both cases, the 
experimental parameters were optimized using the technique described in Section 3.3. 
The thermal property estimates, effective thermal conductivity perpendicular to the fibers 
and effective volumetric heat capacity, obtained for the IM7 graphite fiber - Bismaleimide 
epoxy matrix composite are also discussed. These thermal properties were estimated 
using both the parameter estimation program, MOD BOX, which requires an exact 
temperature solution and the finite element software, EAL, where the temperature solution 
is calculated numerically. In both cases, the properties were estimated using the modified 
Box-Kanemasu method described in Section 3.2. 

The first and second sections of this chapter discuss the optimal experimental 
results obtained for the one-dimensional analysis and the estimated thermal properties 
found utilizing this design, respectively, while the last section discusses the optimal 
experimental results obtained for the two-dimensional configurations. 
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5.1 Results Obtained for the One-Dimensional Analysis (Isotropic Composite 

Material) 

The thermal properties estimated for the one- dimensional analysis include the 
effective thermal conductivity perpendicular to the fiber axis, or the isotropic thermal 
conductivity, and the effective volumetric heat capacity. As mentioned in Section 3.2, 
measured temperatures were required for this estimation procedure; therefore, experiments 
had to be conducted. The next subsections discuss the results obtained for the 
optimization procedure used to determine the optimal experimental parameters utilized in 
these experiments. This includes an analysis of the sensitivity coefficients and the 
determination of the following experimental parameters: the heating time of the uniform 
heat flux, the temperature sensor location, and the total experimental time. 

5.1.1 One-Dimensional Optimal Experimental Design 

The minimization procedure used in this analysis to estimate the effective thermal 
properties is the Box-Kanemasu method. This method requires both measured and 
experimental temperatures. The experiments used to obtain the measured temperatures 
were optimized to provide more accurate property estimates. The optimization technique 
selected in this study, as discussed in Section 3.3, maximizes the determinant of the 
product of the dimensionless sensitivity coefficients and their transpose. Therefore, the 
first step in the optimization procedure is to calculate and analyze the sensitivity 
coefficients for each property. These coefficients are discussed next 
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5. 1 . 1 . 1 Sensitivity Coefficient Analysis 

The first step in the optimization procedure was to calculate and plot the sensitivity 
coefficients. Recall that these coefficients are the derivatives of the temperature with 
respect to the unknown thermal properties and indicate the sensitivity of the temperature 
response due to changes in the parameters. In order for the thermal properties to be 
independently and accurately estimated, these coefficients should be large in magnitude 
(on the order of the temperature rise) and linearly independent. If the sensitivity 
coefficients are small, not enough information is available for the estimation procedure 
and if linear dependence exists between them, the parameters cannot be independently 
estimated. It is important to note, however, that the sensitivity coefficients may be 
linearly dependent over one range but independent over a different range. Knowing these 
ranges will help optimize the physical experiment. 

The dimensionless sensitivity coefficients for the isotropic analysis are given by 
Eqs. (3.52-55). Figures 5.1 and 5.2 show the dimensionless sensitivity coefficients for 
the effective thermal conductivity (X k + ) and effective volumetric heat capacity ( Xj ), 

respectively, at various positions within the composite. In Fig. 5.1, it is seen that after 
a dimensionless time of approximately three, all of the coefficients converge to a constant 
value. This indicates that temperature measurements taken beyond this dimensionless 
time supply little additional information for the estimation of This same result also 
occurs for the effective volumetric heat capacity sensitivity coefficients (Fig. 5.2). Here, 
after a dimensionless time of approximately three, the coefficients converge to zero. To 
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Sensitivity Coefficient, X(kx-eff) 





Sensitivity Coefficient, X(c e jr) 



Figure 5.2. Transient Effective Volumetric Specific Heat Sensitivity Coefficients with the Heat 
Flux, q x ( t), Applied for the Entire Experimental Time for Various x + (x/L x ) 
Locations. 



more accurately estimate C eff> the majority of the temperature measurements should be 
taken over the dimensionless time range of zero to three, where the magnitude of the 
sensitivity coefficients is largest. It is apparent from these figures that the coefficients 
with the largest magnitude occur at a dimensionless x + (defined as x/L x ) location of zero. 
This position corresponds to the heated surface in this analysis. However, it is also 
evident that the X k have a larger magnitude than those for Xq . This indicates that the 

temperature data provides more information about k x _ eff than it does for C eff) and therefore, 
the estimates of k x . eff will be more accurate than those for C eff . 

It should be noted that the sensitivity coefficients in Figs. 5.1 and 5.2 were 
calculated with the heat flux applied for the entire duration of the experiment. Because 
the volumetric heat capacity coefficients approach zero, it suggests that a better scheme 
may consist of applying the heat flux for a finite duration instead of for the entire 
experimental time. This confirms that observing the sensitivity coefficients can give 
insight into the accuracy of the experimental design. 

As mentioned, if the sensitivity coefficients are linearly dependent, the thermal 
properties are correlated and cannot be simultaneously estimated. One way to determine 
if linear dependence exists is to plot the ratio of the sensitivity coefficients, as shown in 
Fig. 5.3 for an x + location of zero. In this figure, the ratio was plotted between the 
dimensionless time range of zero to three. (After this range, the volumetric heat capacity 
sensitivity coefficients converge to zero and the ratio becomes insignificant). If a constant 
curve occurs, the coefficients are linearly dependent. However, as evident from Fig. 5.3, 
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Figure 5.4. Dimensionless Determinant, D + , for Various Dimensionless Heating Times, t h + , as 
a Function of Total Experimental Time, t N + . 



Maximum Determinant, (Dmax) 




(5.1) 



For this study, the thickness of the sample, L x , was 6.78 mm. However, to determine a x , 
k x _ eff and C eff are required, which are the unknown parameters being estimated. Therefore, 
the actual heating time can be estimated by using previous estimates of k x _^ and C eff of 
other similar carbon-epoxy composite materials. The previous estimates used in this study 
were obtained from Scott and Beck (1992a). Using these values, the optimal heating time 
was calculated to be approximately 180 seconds. This value can be updated by 
conducting the experiments, obtaining new estimates for k x _ eff and C eff , recalculating the 
heating time using these new estimates, and repeating the process in an iterative procedure 
until the thermal properties no longer vary. However, as mentioned, a flat peak exists 
between heating times of 2.0 and 2.5. Therefore, the optimal heating time does not have 
to be precise to obtain the most accurate thermal property estimates. 

5. 1.1.3 Optimal Temperature Sensor Location 

Next, the optimal temperature sensor location was determined by plotting the 
determinant as a function of heating time for various sensor locations (Fig. 5.6). It was 
found that the determinant is maximized when the sensor is located at the heated surface 
(x + = 0.0). This result is consistent with the sensitivity coefficients shown in Figs. 5.1 
and 5.2, where the coefficients with the largest magnitude occurred at the heated surface. 
Note that by placing additional thermocouples at other positions within the composite will 
be redundant and will not supply more information for the estimation procedure. 


70 



5 


6 
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Figure 5.6. Determination of the Optimal Temperature Sensor Location, x* (x/L x ). 




5. 1.1.4 Optimal Experimental Time 

The last parameter that was determined was the optimal experimental time. In 
order to see the effect of added data to the value of the determinant, £> + , the determinant 
was calculated from Eqs. (3.42) and (3.43) using the optimal heating time and sensor 
location previously found, but without averaging the integral contained in Eq. (3.43) over 
time. The results are shown in Fig. 5.7; here, it is evident that after a dimensionless time 
of approximately five, the determinant no longer changes significantly. This implies that 
after this dimensionless time, the temperature is reaching its initial state and little 
additional information is being provided for the estimation of the thermal properties. 
Therefore, the experiments can be concluded after a dimensionless time, t N + , of 
approximately five. Note that this is a conservative choice, however, and from Fig. 5.7, 
a smaller value, such as four, could have also been chosen. Again, the actual 
experimental time that t N + represented was found by using the definition of t N + (Eq. 
(3.44c)) and the previous estimates for k x _ eff and C eff (Scott and Beck, 1992a). The 
experimental time calculated for this study was approximately 8 minutes. However, when 
examining the measured temperatures obtained from the experiments, it is seen that its 
initial state (a dimensionless value of zero) is reached after approximately five to six 
minutes, indicating that no new temperature information is being supplied. This 
corresponds to a dimensionless experimental time of three to four, again showing that a 
t N + of five is a conservative value. 
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Modified Determinant, 





5. 1.1.5 Sensitivity Coefficient Using the Optimal Experimental Parameters 

To illustrate the optimal results, the sensitivity coefficients were re-calculated using 
these optimal experimental parameters and are shown in Fig. 5.8 at an x + location of zero. 
As one can see from this figure, more information used to estimate the thermal properties 
is supplied when the heater is applied for the determined optimal heating time rather than 
over the entire experimental time. This occurs because when the heater is applied for the 
entire duration, steady-state values are reached early on and information is no longer 
available for the estimation of Or However, when turning the heater off during the 
experiment, a new transient response is introduced which results in additional temperature 
information for the estimation procedure. 

5.1.2 Estimation of Thermal Properties for Isotropic Materials 

The thermal properties, effective thermal conductivity perpendicular to the fiber 
axis and effective volumetric heat capacity, were estimated for an IM7 graphite fiber - 
Bismaleimide epoxy matrix composite both analytically, using the program MOD BOX, 
and numerically, utilizing the finite element software, EAL. In both cases, the modified 
Box-Kanemasu method was used in the estimation procedure. The properties were 
estimated using both experimental data and numerical and exact solutions so that the two 
could be compared and verification of the accuracy of EAL could be made. The 
estimated thermal properties obtained for the one-dimensional analysis are given in the 
following two subsections. In this analysis, the optimal experimental design previously 
determined was utilized to record the experimental temperatures required. 
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Figure 5.8. Transient Dimensionless Sensitivity Coefficients, X k + and X c + , at an / (x/L x ) 
Location of 0.0 Using the Dimensionless Optimal Heating Time, t h op[ + . 



5. 1.2.1 Estimated Thermal Properties Using an Exact Temperature Solution 

The thermal properties, effective thermal conductivity perpendicular to the fiber 
axis and effective volumetric heat capacity, were estimated from an exact temperature 
model (Eqs. (3.8) and (3.9)) using the sequential, non-linear estimation program 
MODBOX. The estimates that were obtained for three repeated experiments are given in 
Table 5.1, along with their 95% confidence intervals. The confidence intervals for each 
estimated parameter, b { , were approximated by 


b,± 


“(N -p) 


1/2 


‘i-mW - P) 


(5.2) 


where p is the number of parameters estimated, N is the number of data points measured, 
P u is the zth diagonal of the P matrix (Eq. (3.37)) which represents the variance of the 
parameter, S is the sum of the squared residuals, and tj^N-p) is the value of the t 
distribution for (i-ot/2) confidence region and ( N-p ) degrees of freedom (Beck and 
Arnold, 1977). In this study, a considerable number of temperature measurements were 
taken and used in the estimation procedure; therefore, only a slight variance in the 
estimates would be expected. This is in fact the case, as shown by the small confidence 
intervals for each property estimate in Table 5.1. It is also seen that the confidence 
intervals for k x . eff were smaller than those for C eff . This implies that the estimates for k x _ eff 
are more accurate than for This is consistent with the sensitivity coefficients (Figs. 
5.1 and 5.2), where the magnitude of the effective volumetric heat capacity sensitivity 
coefficients is less than that of the effective thermal conductivity. Therefore, the 
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estimation of Or is more sensitive to experimental errors and will not be as accurate as 


estimates for k x . eff . 

The mean value of the thermal property estimates was also calculated, along with 
its 95% confidence interval. In this case, the 95% confidence intervals were obtained 
from 


b, ± L (5.3) 

JN 

where b i and s are the mean and standard deviation of the estimate, respectively, N is 
the number of data points used, and is the value of the t distribution with (N-J) 
degrees of freedom and a/2 confidence region (Walpole and Myers, 1978). As seen in 
Table 5.1 for all three experiments, the property estimates fall within the 95% confidence 
intervals of the mean values. 

To determine how accurately the calculated temperatures matched the measured 
temperatures, the Root Mean Square (RMS) error was also computed where 


E <*; - T ) 2 


RMS = 


Here, T t and Y { are the calculated and measured temperatures, respectively, at the ith time 
step, and N is the total number of temperature measurements. The RMS values were 
calculated two different ways. First, the measured temperatures for each individual 
experiment were compared with calculated values using the thermal properties estimated 
for that experiment; these values are indicated by RMSj in Table 5.1. The RMS values 
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were then determined using the experimental temperatures and calculated temperatures 
determined using the mean thermal property values (also shown in Table 5.1); these 
values are indicated by RMS m . 

To demonstrate the validity of the estimated properties, the calculated temperatures 
obtained using the estimated effective thermal conductivity and effective volumetric heat 
capacity values were compared with the measured temperatures for Experiment 3 in Fig. 
5.9. As one can see, there was very good agreement between the calculated and 
measured temperatures which indicates that the estimated values are reliable. 

The significance of the RMS m values, or the errors resulting from using the mean 


Table 5.1 Estimated effective thermal conductivity, k x ^ and volumetric heat capacity, 
Qip from Experiments 1, 2, and 3, using exact temperature solutions along 
with the Root Mean Square error calculated from individual and mean thermal 
property estimates (RMSj and RMS m ). 



j Exp. 1 

Exp. 2 

Exp. 3 

Mean 1 

K-eff (W/m°C) 

0.519 

±0.002 

0.506 

±0.002 

0.529 

±0.003 

0.518 ±0.028 1 

c ejr (MJ/m 3o C) 

1.423 

1.505 

1.495 

1.474 ±0.111 1 

± 0.013 

± 0.012 

±0.008 

1 

RMSj (°C) 
% Maximum 

0.0526 

0.0815 

0.0652 

mpi iiiii® I 

Temperature Rise 

0.24% 

0.36% 

0.26% 

* :> / % 
S\ •. V "V, ' ' V. 

RMS U (°C) 
% Maximum 

0.0548 

0.0908 

0.0827 

Ss : 

v ■> ■V‘- '■ 

SS -a. -' • \\.\ ■. \-C; 

- ■■ 

Temperature Rise 

0.34% 

0.40% 

0.34% 
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Figure 5.9. Calculated and Mea* 




thermal property estimates, was demonstrated by plotting the temperatures calculated 
using the mean thermal property estimates against each set of experimental data in Fig. 
5.10. From this figure, it is evident that in each case, the calculated temperatures closely 
match the measured temperatures. Furthermore, in comparing Table 5.1 with Fig. 5.10, 
the slight under prediction of temperature in Experiment 2 and the slight over prediction 
of temperature in Experiment 3 can be attributed to the small differences between the 
individual estimated effective thermal conductivity and the mean value, with the thermal 
conductivity estimate for Experiment 2 being slightly under the mean and the value for 
Experiment 3 being slightly over the mean value. However, even with these slight 
variations, the RMS m as a percentage of the maximum temperature rise for each run was 
less than 0.5%, as shown in Table 5.1. This indicated that for all three cases, the mean 
values provided reasonable estimates of the true thermal property values. It also indicates 
that the model used to describe this heat conduction process, as well as the experimental 
design used to obtain the temperature data, are satisfactory. 

The estimate obtained for the effective thermal conductivity was compared with 
results obtained at NASA-LaRC using a cut-bar comparative apparatus (Dynatech, model 
number TCFCM-N4). This device operates by supplying a steady state heat flow in one 
dimension across the composite sample and the same heat flow through a known standard 
material. The temperature difference across the standard material allows for the 
determination of the heat flux, while the temperature difference across the sample gives 
the value of the effective thermal conductivity. Using this apparatus, experiments were 
conducted on three composite samples that were the same type studied in this analysis. 
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Figure 5.10. Temperature (T) Profiles forBoi 
Measurements Using Averagt 




The average result obtained for these three composite samples for the thermal 
conductivity was 0.473 ± 0.038 W/m°C. This is approximately a 9% difference from the 
average k x . eff estimated in this study. However, these experiments were conducted at 
40°C, whereas in this investigation, the experiments were performed at room temperature; 
therefore, an exact comparison cannot be made. 

5. 1.2.2 Estimated Thermal Properties Using a Numerical Temperature Solution 

The thermal properties were also estimated using EAL, where the temperatures 
were calculated numerically. Again, the modified Box-Kanemasu estimation technique 
was employed. The results obtained for the estimated effective thermal conductivity 
perpendicular to the fibers and effective volumetric heat capacity for the three 
experiments are shown in Table 5.2, along with the % difference from the estimates 
obtained using the exact temperature solutions. This % difference is defined as 

where p is the parameter being estimated. As one can see from Table 5.2, the estimates 
found using EAL closely match those obtained using exact temperature models with a 
percent difference of less than 1% occurring. The effective volumetric heat capacity 
estimates had the largest percent differences and resulted because this property is more 
difficult to estimate than the thermal conductivity. As mentioned, this occurs because the 
magnitude of the sensitivity coefficients for C eff are less than those for k Xmeff> causing the 
estimation of C* to be more sensitive to experimental errors and to not be as accurately 
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Table 5.2 Estimated effective thermal conductivity, k x _ e p and volumetric heat capacity, 
from Experiments 1, 2, and 3, using numerical temperature solutions 
(from EAL) along with the % difference from estimates calculated using exact 
temperature models. 


1 Exp. 1 

Exp. 2 

Exp. 3 

Mean j 

K-eff (W/m°C) 

0.518 

0.503 

0.527 

0.516 1 

% Difference 

0.19% 

0.59% 

0.38% 

0.39% 

C#(MJ/m 3o C) 

1.420 

1.495 

1.486 

1.467 

% Difference 

0.21% 

0.66% 

0.60% 

0.49% 


estimated as k x . efr 

From Table 5.2, it can be concluded that the estimated parameters found using the 
finite element software, EAL, are quite accurate and provide reasonable estimates of the 
true thermal property values. 

To verify the accuracy of the temperature solution found using EAL, based on the 
estimated parameters, the temperature profile was plotted along with the temperature 
distribution calculated using an exact analytical solution, as shown in Fig. 11. Here, it 
is seen that the two curves are essentially equal, and therefore, using EAL provides 
reliable temperatures. This is also shown by calculating the RMS value, as given in Eq. 
(5.4), where Y { is the temperature calculated from an exact solution and T { is the 
temperature calculated from EAL. For the experiment shown in Fig. 5.11, the RMS value 
was only 0.27%, again indicating the accuracy of EAL. 
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Figure 5.11. Temperature (T) Profile Using Both an Exact Solution and EAL. 



5. 1.2.3 Sequential Parameter Estimates 

Viewing the sequential parameter estimates can give insight on the validity of the 
mathematical model used to represent the heat conduction process and the resulting 
experimental design. The sequential estimates for the converged values of the thermal 
conductivity and volumetric heat capacity for Experiment 3 are plotted in Fig. 5.12; these 
estimates were obtained using exact mathematical models. From this figure, it is evident 
that each estimate fluctuated greatly towards the beginning of the experiment. This 
occurred because the heat flux had just been activated and not enough temperature 
information was available for the estimation procedure. However, after approximately 
400 seconds which corresponds to a dimensionless experimental time of approximately 
three, the estimates for both the thermal conductivity and volumetric heat capacity are 
constant, indicating that additional data would have provided little additional information 
for the estimation of these parameters. This also indicates that the heat conduction model 
is satisfactory and the optimal experimental time of five is indeed a conservative value, 
as discussed previously. 

5.2 Results Obtained for the Two-Dimensional Analysis (Anisotropic Composite 

Material) 

For the two-dimensional analysis, three properties can be estimated simultaneously: 
effective thermal conductivity perpendicular to the fiber axis (£*_ ejr ), effective thermal 
conductivity parallel to the fiber axis (k y _ e ^ and effective volumetric heat capacity 
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Sequential Estimates 
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Figure 5.12. Sequential Estimates for Effective Thermal Conductivity, k x . eff (W/m°C) and 
Effective Volumetric Heat Capacity, C eff (MJ/m 3o C). 







In this investigation, these thermal properties are not estimated; however, the experimental 
designs used to obtain the temperature data required for the estimation procedure are 
optimized, as in the one-dimensional analysis. Discussed next are the results obtained for 
these optimal experiments. 

5.2.1 Two-Dimensional Optimal Experimental Designs 

Two different two-dimensional experimental configurations were analyzed in this 
study, each containing different boundary conditions. The experimental parameters for 
both configurations were optimized by maximizing the determinant of the product of the 
sensitivity coefficients and their transpose. The maximum determinant values for both 
configurations were then compared to determine which design would be the best choice; 
the configuration with the largest determinant value would give the most accurate 
property estimates. Recall that both configurations had a uniform heat flux imposed over 
a portion of one boundary with the remainder of the boundary insulated. In addition, 
Configuration 1 had known, constant temperatures at the remaining three boundaries, 
while the second configuration had a constant temperature at the boundary opposite to the 
flux boundary and insulated conditions at the remaining two boundaries. (For clarity, 
these configurations are again shown in Figs. 5.13 and 5.14). Therefore, in addition to 
the optimal parameters determined for the one-dimensional case, the optimal position of 
the heat flux was also determined for both configurations. Because of the different 
boundary conditions used, the optimal experimental parameters for each design will not 
be identical. For example, the portion of the boundary that the heat flux should cover 
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Figure 5.13. Experimental Set-up Used for Configuration 1 in the Estimation of the 
Effective Thermal Conductivities in Two Orthogonal Planes. 



Constant 

Temperature 


Figure 5.14. Experimental Set-up Used for Configuration 2 in the Estimation of the 
Effective Thermal Conductivities in Two Orthogonal Planes. 
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will vary between the two configurations, and therefore, must be determined for each 
individual design. The optimal parameter results for both configurations are discussed 
in the following two subsections. These include the optimal heating time, optimal 
temperature sensor location, optimal heat flux position, and optimal experimental time. 
In addition, the sensitivity coefficients for both configurations are analyzed for insight into 
the experimental design and to determine if possible correlation exists between the 
thermal properties. The two configurations are then compared to determine which will 
provide more accurate property estimates, and finally, the last subsection discusses the 
optimal values for various composite dimension (L^) and thermal conductivity (k^) ratios. 

5.2.1. 1 Optimal Experimental Parameters Determined for Configuration 1 

Using Configuration 1, it was desired to select the experimental parameters which 
maximize the sensitivity of the temperature with respect to all of the unknown thermal 
properties. The same technique was used as in the one-dimensional optimization 
procedure, only now, the required maximum temperature value (T mai + ) was redefined as 
the temperature attained at steady state. This T BUU + was used because it represents the 
actual maximum temperature that could be reached for this particular design. 

To perform the optimization procedure, the temperature solutions and sensitivity 
coefficients require predetermined values for L xy (LJLy) and (k y _ e£ /k x . eff ). The value 
chosen for L xy in this analysis corresponds to the size (0.49 cm x 10.16 cm) of an existing 
composite sample that can be used in the experiments to determine the temperatures 
needed in the estimation procedure, while the value for was taken from previously 
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measured effective thermal conductivities parallel and perpendicular to the fiber axis of 
similar carbon-epoxy composite samples (Loh and Beck, 1991). The specific values used 
were L xy - 0.048 and - 1. However, to allow the optimization procedure to be 
applicable to other composite dimensions or effective thermal conductivity ratios, optimal 
experimental parameters were also calculated for all possible combinations of L xy s equal 
to 0.5 and 1.0, and x^’s equal to 1 and 1/7. This results in a total of nine combinations. 
The results for the combination discussed above (L xy = 0.048 and k x> , = 7) will be 
examined the most thoroughly, however, since these are the actual conditions of an 
existing composite sample that can later be utilized in the experimental designs to 
estimate the thermal properties. 

In performing the optimization procedure for L xy = 0.048 and K xy = 7, five 
parameters were optimized: the dimensionless portion of the boundary that the heat flux 
is applied, the dimensionless location (x s + 9 y s + ) of the temperature sensor, the 

dimensionless heating time, t h + , and the dimensionless experimental time, t N + . Note that 
the optimal experimental time is not as important as the other four parameters. Therefore, 
the optimal procedure used to determine x s + , y/, L p l + , and t h + did not take into account 
the optimization of t N + (this value was determined last). 

The most accurate way to determine the optimal value for each parameter is to 
differentiate the determinant given in Eq. (3.45) with respect to each of the experimental 
parameters and set it equal to zero, resulting in four equations and four unknowns. These 
equations can then be solved simultaneously, allowing for the desired optimal parameters 
to be determined. However, this method is not practical due to the complexity of the 
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equations involved. Therefore, an iterative scheme was developed where a program was 
used to vary each of the four parameters individually (excluding t N + ). This iterative 
procedure consists of two phases; the first phase includes determining the general range 
of the optimal values, while the second phase narrows this range to determine the optimal 
experimental parameters more precisely. In phase one, the following procedure was used: 

1) Fix x s *, y*, and L p / to their starting values (0.0, 0.0, and 0.1, respectively). 

2) Vary t h + from 0.05 to 5.0 by 0.05. For each t h + , calculate the determinant, D + , as 
a function of time, t + . 

3) Determine the maximum determinant value, D max + , for each of the determinant 
curves generated in step 2 for each t h + . 

4) Compare the maximum determinant values found for each t h + and record the one 
with the largest magnitude, along with its corresponding heating time, t h + . 

5) Holding y s + and L p / constant at their original values, vary x + and repeat steps 2 
through 4. Note, x s + was varied from 0.0 to 1.0 by 0.1 increments. 

6) After the x s + loop is completed, change y s + to its new value (increment the previous 
value by 0.1) holding 1 fixed and again repeat steps 2 through 5. Note, y* was 
varied from 0.0 to 1.0 by 0.1 increments. 

7) Finally, change L pI + to its new value (increment the previous value by 0.1) and 
repeat steps 2 through 6 in the designated order. 

This procedure then provides a maximum determinant value for all combinations of x/, 
y s + and Lp/, with the corresponding t h + . From this data, the general region of the actual 
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maximum determinant can be determined. Phase two then involves refining the grid sizes 
for the parameters in this region to determine the location more precisely. Since 
there is more than one parameter that can vary, the procedure is more complex than the 
one-dimensional analysis, and must be iteratively updated. 

5.2. 1.1.1 Optimal Temperature Sensor Location on the x + Axis 

The first optimal parameter determined was the temperature sensor location along 
the x + direction. Recall from the experimental configuration (Fig. 5.13) that this is the 
direction parallel to the heat flow. From the one-dimensional analysis, it was determined 
that the optimal location to place the sensor was at the heated surface. Therefore, the 
same result would be expected for the two-dimensional analysis. This was in fact the 
case, with the maximum determinant always occurring at a x s + location equal to zero (or 
at the heat flux boundary) for all combinations of Lp/ and y+. This result is reasonable 
because the maximum determinant occurs when the sensitivity coefficients are the largest, 
or when the greatest temperature variation occurs. Since the temperature of the composite 
is initially at a dimensionless value of zero, then at the boundary where the uniform heat 
flux is applied would be the location where the largest temperature gradient in the x + 
direction would occur. 

5.2.1. 1.2 Optimal Temperature Sensor Location on the y + Axis 

The next parameter chosen to optimize was the sensor location along the y + axis. 
After calculating the maximum determinant value for each x s + , y s + , and 1^/ combination 
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and determining that the optimal sensor location along the x + axis was at the heated 
surface, it was found that the true maximum determinant (the largest value of all of the 
maximum determinants for each combination) was in the general region of y/=0.1, 
Lp/- 1-0, and t h *-\35. However, it should be noted that for all L pl + locations, the 
maximum determinant always occurred at y* equal to 0.1 with a t h + of approximately 
1.35. Using the optimal values of L p /- 1.0 and 4 + =1.35, the grid size for y* was refined 
around 0.13, using a range from 0.05 to 2.0, to determine the optimal y* location more 
precisely. Using this refined range, the maximum determinant occurred at a new y* 
location of 0.13, as shown in Fig. 5.15, where the maximum determinant values for 
various y s + locations are plotted (again, using Lp/= 1.0 and 4 + =1.35). It should be noted 
that when the heat flux is applied across the entire boundary (L p /= 1.0), the problem 
becomes symmetric. Therefore, a y* of 0.87 would also be an optimal location, resulting 
in the same maximum determinant value as for y s + equal to 0.13. 

5.2. 1.1.3 Optimal Heating Time 

The next parameter determined was the optimal heating time. Using the optimal 
location for y s + found above of 0.13 and the corresponding L pl + of 1.0, the dimensionless 
heating times were varied around the previously determined optimal value of t h + =135 
(ranging from 0.05 to 2.0). The maximum determinant then occurred at a new t h + of 1.4, 
as seen in Fig. 5.16, where the maximum determinants are plotted for various t h + values. 
As mentioned, since more than one parameter can vary, determining the optimal 
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Figure 5.15. Maximum Determinant, D max + , for Various Locations Along the y + (ylL y ) Axis 
Calculated Using the Optimal Values of x s + =0.0, L p] + - 1.0, and t h + =l.35. 
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Figure 5.16. Maximum Determinant, D max + , for Various Heating Times, t h + . Calculated Using 
the Optimal Values of x/=0.0, L p /= 1.0, and y/=0.13. 



parameters for the actual maximum determinant becomes an iterative process. Therefore, 
since a new t h + was determined, the y + values were again varied over the same range 
(0.05 to 2.0) using the new heating time of 1.4 (and the previously calculated optimal 
value for Lp / of 1.0) to see if its optimal value changed. However, as seen in Fig. 5.17, 
changing the heating time from 1.35 to 1.4 did not alter the optimal y s + value of 0.13. 

5.2. 1.1. 4 Optimal Heat Flux Location, L,/ 

Using the optimal parameters determined of */=0.0, y/=0.13, and t h + = 1.4, the 
position of the heat flux, L pl + y was then varied from 0.6 to 1.0 to see if the previous 
optimal location of Z^ /=1.0 changed when using these new y s + and t h + values. This 
result is shown in Fig. 5.18. As seen in this figure, the maximum determinant occurred 
at a Lp/ location of 1.0, as obtained previously. This indicates that the optimal design 
for Configuration 1 consists of having the heat flux applied over the entire boundary. 
However, it is evident that the curve in Fig. 5.18 is rather flat when the heater is applied 
over 70 to 100% (0.70 to 1.0) of the boundary, and therefore, any value in this range 
could be used to obtain the same accuracy in the property estimates. 

5 .2. 1 . 1 .5 Optimal Experimental Time 

Finally, the last parameter determined was the optimal dimensionless experimental 
time, t N + . This was calculated using the same procedure as for the one-dimensional 
analysis, where the dimensionless determinant, Z) + , was calculated from Eqs. (3.42) and 
(3.43) using the optimal parameters determined for x+, y/ L pl + , and t h + , but without 
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Figure 5.17. Maximum Determinant, D max + , for Various Locations Along the y + (y!L y ) Axis 
Calculated Using the Optimal Values of x/=0.0, L pl + = 1.0, and ^ + =1.35 (Second 
Iteration). 



Maximum Determinant, (Dmax) + x 10 




averaging the integral contained in Eq. (3.43) over time. The results are shown in Fig. 
5.19; here, it is evident that after a dimensionless time of approximately 4.0, the 
determinant no longer changes significantly. This implies that after this dimensionless 
time, the temperature is returning to its initial state (a dimensionless value of zero) and 
little additional or no information is being provided for the estimation of the thermal 
properties. Therefore, the experiments can be concluded after a dimensionless time, t N + , 
of approximately 4.0. Again, however, as in the one-dimensional case, this is a 
conservative choice, and a smaller value, such as 3.5, could have also been chosen. 

5.2. 1.1.6 Verification of the Optimal Temperature Sensor Location of the x + Axis 

To verify the optimal location of the temperature sensor along the x + direction, for 
which a value of zero was determined, the dimensionless determinant was calculated 
using the optimal values for y s + , L p } + , and t h + for various x a + locations. The results are 
plotted in Fig. 5.20 against dimensionless time. As seen from this figure, the maximum 
determinant occurred when the sensor was at the heated surface (x s + = 0.0), confirming 
the optimal result obtained for the x* location. 

5.2. 1.1.7 Maximum Determinant Using the Optimal Experimental Parameters 

In summary, the above optimization procedure resulted in the following optimal 
experimental parameters for Configuration 1: x/=0.0, y/=0.13, L pI + =l.0, and t h + = 1.4. 
Using these optimal experimental parameters, the dimensionless determinant, D + , was 
plotted versus dimensionless time, where a maximum of 5.36 x 10 7 occurred, as shown 
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in Fig. 5.21. The reason why these determinant values are less than those obtained for 
the one-dimensional case is because D + is now a 3 x 3 determinant, as given in Eq. 
(3.45). Since the sensitivity coefficients are of the same order of magnitude for both the 
one-dimensional and two-dimensional cases, then multiplying three coefficients together, 
as required in the 2-D determinant, will result in smaller maximum determinant values 
than multiplying only two coefficient values, as in the 3-D determinant 

5.2.1. 1.8 Temperature Distributions for Configuration 1 

Using the optimal values for x s + and L pJ + , temperature was plotted for various y* 
locations for four different dimensionless times: early (0.1), intermediate (two at 0.5 and 
1.4), and steady state (4.0) (Fig. 5.22). Note that these temperature distributions were 
calculated with the heat flux applied for the entire experimental time, t N + . The desired 
optimal values occur when the determinant is a maximum, or when the sensitivity 
coefficients are the most sensitive to temperature changes. This typically occurs when 
the temperature gradient is large. As seen from this figure, at and near the optimal y* 
location of 0.13, the temperature gradient is steep with respect to y s + , and therefore, the 
sensitivity coefficients for k y . eff are expected to be large in magnitude. This steep gradient 
occurs because the composite sample is heated, however, the temperatures at the 
boundaries are held constant, resulting in a large temperature variation. However, it is 
also seen that the optimal y* location is not at the steepest temperature gradient. This 
results because the maximum determinant occurs when the product of the sensitivity 
coefficients for all three parameters is the largest in magnitude, which is not necessarily 
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Figure 5.21. Dimensionless Determinant, D + , Calculated Using the Optimal Experimental 
Parameters of * 5 + =0.0, j s + =0.13, L p /= 1.0, and t h + =lA. 
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Figure 5.22 Temperature ( T) Distribution for Various y* Locations Calculated Using Four 
Different Experimental Times. 



at the largest temperature gradient in the y + direction. 

The temperature distribution was also calculated as a function of time using the 
optimal experimental parameters determined above (see Fig. 5.23). As shown in Fig. 
5.23, the heat flux is terminated as the temperature approaches steady state (t h + = 1.4). 
This is consistent with Fig. 5.22, where applying the heat flux for the dimensionless time 
of 1.4 results in temperatures close to the steady state temperatures attained at C=4.0. 

5.2. 1.1.9 Sensitivity Coefficients Calculated Using the Optimal Experimental Parameters 

Using the optimal experimental parameters determined, the dimensionless 
sensitivity coefficients for the three effective thermal properties, k re p and C eff) were 
calculated and plotted as a function of dimensionless time, as shown in Fig. 5.24. Here, 
it is seen that the sensitivity coefficients for k x _ eff and C^are relatively large in magnitude, 
being on the same order as the temperature rise, with the k x . eff coefficients being the 
largest. The sensitivity coefficients for k y _ eff have the smallest magnitude of all three. It 
is also seen that after a dimensionless time of approximately four, all of the coefficients 
converge to zero, indicating that temperature measurements taken beyond this time supply 
little additional information for the estimation procedure. This result is consistent with 
the temperature distribution (Fig. 5.23), where its initial state was attained after this 
dimensionless time. Therefore, no new temperature information is being provided and 
the estimation procedure is complete. This result is also consistent with the determined 
optimal experimental time, where after a t N + of 4.0, the determinant no longer varied. 

Since the sensitivity coefficients for k x _ eff have a larger magnitude than the 
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Figure 5.23. Temperature Distribution for Configuration 1 Calculated Using the Optimal 
Experimental Parameters of x s + =0.0, y/=0.13, L p /= 1.0, and t h + = 1.4. 
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Figure 5.24. Dimensionless Sensitivity Coefficients, XC ,and,X c + , Calculated 
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Using the Optimal Experimental Parameters of x/=0.0, y, + =0. 13, L p , + =1.0, and 
V= 1.4. 



coefficients for k y . eff and C eff> it implies that the temperature data are supplying more 
information for the estimation of k x . eff than k y . eff and C eff . As a result, the estimated values 
obtained for k x _ eff can be regarded as the most accurate of the three parameters, resulting 
in the smallest confidence intervals. 

As mentioned in the one-dimensional case, it is important to plot the sensitivity 
coefficients to see if they are correlated. If correlation occurs, the thermal properties 
cannot be estimated independently. From Fig. 5.24, it is evident that the Q sensitivity 
coefficient, which changes from negative to positive values, is not correlated with either 
the k x . eff or k y . eff coefficients, which are always negative. However, this observation is not 
as apparent between the k x _ eff and k reff sensitivity coefficients. If they are correlated, then 
only the ratio, k^,, can be estimated. Therefore, to test for possible correlation, the ratio 
of XC /Xl was calculated and plotted as a function of dimensionless time (Fig. 5.25). 

rf ff *-*ff 

If a straight line occurs, the two parameters are correlated. However, as evident from Fig. 
5.25, the line is far from linear, and therefore, k^ and k y _ eff can be estimated 
simultaneously. 

5.2. 1.2 Optimal Experimental Parameters Determined for Configuration 2 

Recall that Configuration 2 consisted of a uniform heat flux imposed over a 
portion of one boundary, with the remainder of the boundary insulated. The boundary 
opposite to the flux boundary was maintained at a known constant temperature, and the 
remaining two boundaries were insulated (Fig. 5.14). Again, for this configuration, it was 
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desired to select the experimental parameters which maximize the sensitivity of the 
temperature with respect to all of the unknown thermal properties. Since the same 
composite samples will be used in the experiments for both Configurations 1 and 2, the 
result for a L ^ of 0.048 and a of 7 (Loh and Beck, 1991) will again be the most 
thoroughly analyzed. However, as in the Configuration 1 case, all possible combinations 
of k ^ equal to 1 and 1/7, and equal to 0.048, 0.5, and 1.0 will also be performed. The 
same optimization procedure, as discussed in Section 5.2.1. 1 for Configuration 1, was 
used and similar experimental parameters were optimized (x s + , y+ 9 L p2 + , t h + , and t N + ). 
These experimental parameter results are discussed next 

5.2. 1.2.1 Optimal Temperature Sensor Location on the x + Axis 

The first optimal experimental parameter determined was the temperature sensor 
location along the x + axis. Again, as in Configuration 1, the maximum determinant 
always occurred at a x* location equal to zero (or at the heated surface) for all 
combinations of y s + and L p,2 + ‘ This occurs for the same reason as discussed in Section 
5.2.1. 1.1, where the largest temperature gradient in the x + direction occurs at the heated 
surface. 

5.2. 1.2.2 Optimal Temperature Sensor Location on the y + Axis and Heat Flux Position 
The next parameters that were optimized were the sensor location along the y + axis 

and the position of the applied heat flux, L p2 + . After calculating the maximum 
determinant for each x* 9 y s + , and L p / combinations, the true maximum determinant (the 
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largest value of all of the maximum determinants for each combination) was in the 
general region of L p2 + = 0.9, y/=0.8, and t h + =1.55. To determine the precise location of 
the actual maximum determinant, both the Lp/ and y s + grid sizes were refined around the 
previously obtained values of 0.9 and 0.8 respectively, (L p 2 + was varied from 0.84 to 0.92 
and y s + was varied from 0.6 to 0.9) while holding t h + constant at 1.55. For each L p2 
value, the maximum determinant was plotted as a function of the y s + location in Fig. 5.26. 
As seen in this figure, the actual maximum determinant occurs at a L p2 + of 0.89 and a 
corresponding y s + location of 0.77. However, it is seen that the curve is fairly flat when 
Lp 2 is located between 0.88 to 0.9; therefore, any value within this range could be used 
for Lp 2 to improve the accuracy of the property estimates. Note that L p2 is different than 
L p ] + , as expected. If L p2 had equalled L pJ + the heat flux would be applied over the 
entire boundary. Due to the insulated boundary conditions on the sides used in this 
configuration, the problem would reduce to one-dimensional heat conduction and k y _ eff 
could no longer be estimated. 

5.2. 1.2.3 Optimal Heating Time 

The next parameter optimized was the heating time. Setting L p 2 and y* equal to 
their optimal values calculated in the above section, (0.89 and 0.77, respectively) the 
heating time was varied around its previously determined value of t h + = 1.55 (from 1.45 to 
1.65). At each heating time, the dimensionless determinant was calculated as a function 
of dimensionless time. A few of these determinant curves are shown in Fig. 5.27 for 
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Figure 5.26. Maximum Determinant, D max + , for Various y s + (y/L y ) Locations and Heater 
Locations, L p>2 + (Lp JL y ), Calculated at x s + =0.0. 
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Figure 5.27. Maximum Determinant, D max + , for Various Heating Times, t h + . Calculated Using 
the Optimal Experimental Parameters of x/=0.0, y s + =0Jl, and L p 2 + =0.89. 



various heating times. The maximum determinant value of all of these curves occurred 
when the heating time was equal to 1.55. This is the same result obtained previously, and 
therefore, the optimal values calculated for 2 + and y s + , which were found using a t h + of 
1.55, do not have to be iteratively updated. Note, however, that the maximum 
determinants are practically equal for all heating times between 1.45 and 1.65. Therefore, 
using a t h + of 1.55 does not have to be precise to provide the most accurate property 
estimates. 

5.2. 1.2.4 Optimal Experimental Time 

Finally, the last parameter determined was the optimal dimensionless experimental 
time, t N + . This was calculated using the same procedure as in Configuration 1 (Section 
5.2. 1.1. 5) with the modified dimensionless determinant results shown in Fig. 5.28. Here, 
it is evident that again, as in Configuration 1, after a dimensionless time of approximately 
4.0, the determinant no longer changes. This implies that after this dimensionless time, 
the temperature is returning to its initial state (a dimensionless value of zero) and little 
additional information is being provided for the estimation of the thermal properties. 
Therefore, the experiments can be concluded after a t N + of 4.0. Again, however, as was 
the case for Configuration 1, this is a conservative choice and a smaller value, such as 
3.5, could have also been chosen. 

5.2. 1.2.5 Verification of the Optimal Temperature Sensor Location on the x + Axis 

To verify the optimal location of the temperature sensor along the x + direction, for 
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which a value of zero was determined, the dimensionless determinant was calculated 
using the optimal values for y s + , L p 2 + , and t h + for various x s + locations. The results are 
shown in Fig. 5.29 against dimensionless time. As seen from this figure, the maximum 
determinant occurred when the sensor was at the heated surface (x/=0.0), confirming the 
optimal result obtained for the x s + location. 

5.2. 1.2.6 Maximum Determinant Using the Optimal Experimental Parameters 

In summary, the above optimization procedure resulted in the following optimal 
experimental parameters for Configuration 2: x/=0.0, y 5 + = 0.77, L p>2 + = 0.89 and t h + = 1.55. 
Using these optimal parameters, the dimensionless determinant, D + , was plotted as a 
function of dimensionless time, where a maximum of 4.29 x IQ 7 occurred, as shown in 
Fig. 5.30. Again, the reason why these determinant values are less than those obtained 
for the one-dimensional case is the same as discussed in Section 5.2. 1.1.7. 

5.2. 1.2.7 Temperature Distributions for Configuration 2 

Using the optimal values for x s + and L p /» the temperature was plotted for various 
y s + locations for four different dimensionless times; initial (0.1), two intermediate (0.5 and 
1.55), and steady state (4.0) (Fig. 5.31). Note that these temperature distributions were 
calculated with the heat flux applied for the entire experimental time, t N + . The desired 
optimal values occur when the determinant is a maximum, or when the sensitivity 
coefficients are the most sensitive to temperature changes. As mentioned previously, in 
the case of thermal conductivity, this can occur when the temperature gradient is large. 
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Figure 5.30. Dimensionless Determinant, D + , Calculated Using the Optimal Experimental 
Parameters of x/=0.0, y/=0.77, L p2 + = 0.89, and f A + =1.55. 
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Figure 5.31 Temperature (T) Distribution for Various y+ Locations Calculated Using Four 
Different Experimental Times. 



As seen from this figure, at an optimal y s + location of 0.77, the temperature gradient is 
steep, and therefore, it is expected that the sensitivity coefficients for k y _ eff are large in 
magnitude. However, the optimal y s + location for this design again does not occur at the 
steepest temperature gradient and results for the same reasons given in the Configuration 
1 analysis. 

The temperature was also calculated as a function of dimensionless time using the 
optimal experimental parameters determined previously (Fig. 5.32). As seen in Fig. 5.32, 
the temperature distribution behaves the same way as for Configuration 1, where the heat 
flux is terminated as the temperature approaches steady state (t h + = 1.55). This is again 
consistent with Fig. 5.31, where applying the heat flux for the dimensionless time of 1.55 
results in temperatures close to the steady state temperatures attained at t + = 4.0. 

5.2. 1.2.8 Sensitivity Coefficients Using the Optimal Experimental Parameters 

Using the optimal experimental parameters determined, the dimensionless 
sensitivity coefficients for the three effective thermal properties, k x ,^, k y . eff , and were 
calculated and plotted as a function of dimensionless time in Fig. 5.33. Here, it is seen 
that the sensitivity coefficients for k x . eff and C eff are relatively large in magnitude, while 
for k y _ eff> the coefficients are much smaller (of the order 0.1). It is also seen that after a 
dimensionless time of approximately four, all of the coefficients converge to zero, 
indicating that temperature measurements taken beyond this time supply little additional 
information for the estimation procedure. This result is consistent with both the 
temperature distribution in Fig. 5.32, where its initial state was attained after this 


120 








Dimensionless Sensitivity Coefficients 



Dimensionless Experimental Time, ( hj + 


Figure 5.33. Dimensionless Sensitivity Coefficients, X^X^and,X Cff , Calculated 
Using the Optimal Experimental Parameters of x 5 + =0.0, y s + =0.77, Z^ 2 + =0.89, and 
4 + =1.55. 



dimensionless time, and with the determined optimal experimental time, where after a t N + 
of 4.0, the determinant no longer varied (Fig. 5.28). 

Again, the sensitivity coefficients should be analyzed to see if they are correlated. 
If correlation occurs, the thermal properties cannot be estimated independently. From Fig. 
3.33, linear independence is again evident between whose coefficient changes from 
negative to positive values, and both k x _ pff and k y . eff) where the coefficients are always 
negative. However, as with Configuration 1, it was desired to determine if the sensitivity 
coefficients for k x . eff and k y . eff are correlated. If correlation occurs, then only the ratio, K^, 
can be estimated. To test for possible correlation, the ratio of X k JX k ^ was again 

calculated, with the results shown in Fig. 5.34. If a straight line occurs, correlation exists. 
From this figure, however, it is evident that a linear line does not occur, and therefore, 
the coefficients are linearly independent and k x _ eff and k y . eff can be estimated 
simultaneously. 

5.2. 1.3 Comparison of Configurations 1 and 2 

After calculating the optimal experimental parameters for both Configurations 1 
and 2, the configurations were compared to determine which one would provide the most 
accurate thermal property estimates. This comparison can be made by determining which 
configuration has the largest maximum determinant. Figure 5.35 shows the dimensionless 
determinants as a function of dimensionless experimental time calculated using the 
optimal experimental parameters found for each configuration. From this figure, it is 
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Figure 5.35. Comparison of the Dimensionless Determinant, D + , for Configurations 1 and 2. 



evident that the design with constant temperatures on the two sides (Configuration 1) will 
provide more accurate estimates than the design with insulated sides (Configuration 2). 
However, this conclusion requires further analysis by viewing the sensitivity coefficients 
for each configuration, as shown in Fig. 5.36. Here, it is seen that the sensitivity 
coefficients for k x _ eff and C eff are 20% and 9% larger in magnitude, respectively, for 
Configuration 2 than Configuration 1. However, the sensitivity coefficient for k y _ eff is 
136% larger in magnitude for Configuration 1 than Configuration 2. (Note, these percents 
were calculated using the largest point on each of the curves). This difference in k y . eff is 
much more substantial than that for and C eff . From viewing these sensitivity 
coefficients, it can be concluded that when estimating all three thermal properties 
simultaneously, Configuration 1 should be utilized, since it will provide approximately 
the same amount of information for k x _ eff and C eff that Configuration 2 would provide, but 
considerably more information for the estimation of k y _ efr This result seems reasonable 
since a greater temperature variation would occur in the y + direction (the same direction 
as k reJ f) when the walls are maintained at a constant temperature of zero rather than 
insulated, where the wall temperatures are allowed to rise (only the gradient at the wall 
is required to remain equal to zero). This result in consistent with comparing the 
maximum determinant values between the two configurations, as discussed previously. 

5.2. 1.4 Other Optimized Parameters 

For both configurations used in the two-dimensional analysis, the optimal 
parameters were determined using a L xy of 0.048 and a K xy of 7. A L xy ratio of this 
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Figure 5.36. Comparison of the Dimensionless Sensitivity Coefficients, X k , X k , and, X^, for 
Configurations 1 and 2. 




magnitude is typical of the sizes of composite samples used in experimental designs. The 
samples used by Loh and Beck (1991) to determine the effective thermal conductivities 
parallel (k y _ eff ) and perpendicular (k x _ eff ) to the fibers, from which they determined a thermal 
conductivity ratio of 7, were also of this magnitude. Since this study uses similar carbon- 
epoxy composite materials, a of 7 was also used in this investigation. However, to 
demonstrate how this optimization analysis could be extended to other composite 
dimensions or effective thermal conductivity ratios, different values for and L xy were 
also used in the optimization procedure. These combinations include L xy equal to 0.048, 
0.5, and 1.0, and K xy equal to 7, 1, and 1/7. The results for all combinations are discussed 
in the following subsections. 

5.2. 1.4.1 Various L ^ and Combinations Used for Configuration 1 

The first combination investigated was = 0.5, and K xy = 7. This L xy results in 
the composite thickness in the direction of heat transfer (the x + direction in this analysis) 
being ten times greater than when L xy equalled 0.048. Using the same procedure 
discussed in Section 5.2. 1.1, the general region of the maximum determinant was 
calculated. Using the optimal experimental parameters found for this region, the 
sensitivity coefficients were calculated and plotted, as shown in Fig. 5.37. From this 
figure, it is seen that the coefficients reach steady state very quickly. This occurs because 
of the significant increase in the thickness in the x + direction, creating more material to 
absorb the heat produced from the applied heat flux. Therefore, it can be concluded that 
using this L xy and combination provides inadequate information for the estimation 
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procedure and is not recommended as an experimental design. Based on this result, it is 
evident that raising L xy to 1.0, where the thickness in the x + direction is additionally 
increased, will provide similar results, and therefore, should also not be used as an 
experimental design. When comparing these two results to the case previously analyzed 
in Section 5.2.1. 1.9 (. L ljcy = 0.048, K ^ = 7), it is seen that when using a thin sample for this 
thermal conductivity ratio, more information is available for the estimation of the thermal 
properties. This is shown by the larger sensitivity coefficients that result for all three 
parameters (Fig. 5.24). 

Next, to determine the effects of different effective thermal conductivity ratios, 
was decreased to 1, and again, combinations for L xy equal to 0.048, 0.5, and 1.0 were 
analyzed. Note that a of 1 implies that the resistance to heat flow is equal in both the 
x + and y + directions (due to equal effective thermal conductivities). 

For all three L xy -K xy combinations, the sensitivity coefficients were plotted using 
experimental parameters around the maximum determinant region. These results are 
shown in Figs. 5.38, 5.39, and 5.40 for L xy values of 0.048, 0.5, and 1.0, respectively. 
From Fig. 5.38, it is seen that when L xy equals 0.048, sufficient information is provided 
for the estimation of k x . eff and C eff> where the sensitivity coefficients are large in 
magnitude. However, the coefficients for k y _ eff are quite small, remaining practically zero 
for the entire experimental time. This implies that the estimation of k y . eff will be difficult, 
and most likely, inaccurate. 

For L xy equal to 1.0 (Fig. 5.40), it is again seen that the sensitivity coefficients 
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reach steady state fairly rapidly, and therefore, little information is being supplied for the 
estimation of the thermal properties. 

The sensitivity coefficients calculated using a L xy of 0.5 (Fig. 5.39), however, are 
all relatively large in magnitude. This implies that when using this L xy ratio, difficulty 
in the estimation of the thermal parameters will not be encountered. 

In conclusion, when a composite has equal effective thermal conductivities parallel 
and perpendicular to the fibers, the optimal L xy ratio is not for either a real small or real 
large thickness in the direction of heat flow, but instead, falls somewhere in between. An 
L xy of 0.5 may perhaps be the optimal ratio; however, this conclusion would require 
further analysis. 

The last combination of L xy s and K* y for Configuration 1 were again, L xy equal to 
0.048, 0.5, and 1.0, with equal to 1/7. Now, the effective thermal conductivity 
perpendicular to the fibers (in the direction of the heat flow) is 7 times larger than the 
thermal conductivity parallel to the fibers. The sensitivity coefficients for all three 
combinations, L xy equal to 0.048, 0.5, and 1.0, were calculated using the optimal 
experimental parameters determined around the maximum determinant region. These are 
shown in Figs. 5.41, 5.42, and 5.43, respectively. From Fig. 5.41, where L xy = 0.048, it 
is seen that k yeff cannot be estimated since the sensitivity coefficient is zero. This occurs 
because the larger thermal conductivity (k x . eS ) is parallel to the heat flow (in the x + 
direction). Since the sample is so thin in this direction, the majority of heat is conducted 
along this path, causing very small temperature variations to occur perpendicular to the 
heat flow, along the y + axis. 
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From Fig. 5.42, is seen that increasing L xy to 0.5 provides better results for the 
estimation of k y . eff) where the sensitivity coefficient is larger in magnitude. This occurs 
because the thickness in the direction of heat flow, x + y increases. Therefore, the 
conduction process is slowed, allowing more heat to be dissipated in the y + direction. 

This result is even more significant when increasing L xy to 1.0 (Fig. 5.43), where 
the largest magnitude for the k y _ eff sensitivity coefficient out of all three combinations 
occurs. 

It can be concluded from these results that when the effective thermal conductivity 
in the direction of heat flow is much larger than the thermal conductivity in the direction 
perpendicular, a better experimental design would consist of a larger thickness in the jc + 
direction, allowing more heat to be dissipated in the direction perpendicular to the heat 
flow. 

5.2. 1.4.2 Various L xy and Combinations Used for Configuration 2 

Again, the first combination investigated was L xy = 0.5 and = 7. The sensitivity 
coefficients calculated using the optimal experimental parameters in the region of the 
maximum determinant are shown in Fig. 5.44. Here, it is seen that not much information 
is being supplied for the estimation of k y _ e p where the sensitivity coefficient is small. 
This occurs because of the increased thickness in the x + direction. Raising L xy to 1.0 will 
further increase this thickness, and therefore, similar results are expected. Therefore, it 
can be concluded that when the thermal conductivity in the direction of heat flow is much 
smaller than the thermal conductivity perpendicular to the direction of heat flow, a 
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composite sample with a small thickness should be used for the best results. This result 
is consistent with the case previously analyzed (L^ = 0.048, = 7), where the sensitivity 

coefficients for all three parameters are slightly larger than when = 0.5. 

Next, was decreased to 1 and similar combinations for were analyzed. 
Again, note that this implies that the resistance to heat flow is equal in both the x + and 
y + directions. The sensitivity coefficients calculated using the experimental parameters 
around the maximum determinant region for all three L^-k x> , combinations are shown in 
Figs. 5.45, 5.46, and 5.47. In Fig. 5.45 (L^ = 0.048), it is seen that having a small 
thickness in the x + direction creates difficulty in the estimation of k y _ eff> where the 
sensitivity coefficient is essentially zero for the entire experimental time. This result is 
consistent with that obtained for Configuration 1 at a similar ratio. However, when 
L xy is increased to 0.5 or 1.0 (Figs. 5.46 and 5.47, respectively), the k y _ eff sensitivity 
coefficients are of approximately the same magnitude, only different in sign. However, 
this magnitude is still small and therefore, neither an of 0.5 or 1.0 is the optimal 
value. The optimal L xy may lie between these values, but the exact determination would 
require further analysis. 

The last combinations of L xy s and for Configuration 2 were again, equal 
to 0.048, 0.5, and 1.0, with K xy equal to 1/7. Now, the effective thermal conductivity 
perpendicular to the fibers (in the direction of heat flow), k x . eff , is 7 times larger than the 
thermal conductivity parallel to the fibers, k y _ eff . The sensitivity coefficients for all three 
combinations, equal to 0.048, 0.5, and 1.0, were calculated using the optimal 
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Figure 5.46. Sensitivity Coefficients for Configuration 2 Using a L xy ( LJL y ) of 0.5 and 
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experimental parameters determined around the maximum determinant region. These 
results are shown in Figs. 5.48, 5.49, and 5.50, respectively. From each of these figures, 
the same results are seen as for Configuration 1, where the best experimental design 
consisted of a larger thickness in the x + direction (a L xy of 1.0), allowing more heat to be 
dissipated in the direction perpendicular to the applied heat flux. This is evident by the 
larger sensitivity coefficient for k y€ff when L xy = 1.0 (Fig. 5.50) than when L xy equals either 
0.5 (Fig. 5.49), where the magnitude of the coefficient is 0.1, or 0.048 (Fig. 5.48), where 
the coefficients are essentially zero. Recall that this zero coefficient results for L rv = 
0.048 because the sample is very thin in the direction of heat transfer, x + (and the 
direction of the larger effective thermal conductivity). Therefore, the majority of the heat 
is conducted along the x + axis, causing very small temperature variations to occur along 
the y + axis, and as a result, the estimation of k y _ eff becomes difficult. 
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Figure 5.48. Sensitivity Coefficients for Configuration 2 Using a L xy ( LJL y ) of 0.048 and 
cyy of i/7. 
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Chapter 6 


Conclusions and Summary 

The primary objectives of this study were to develop optimal experimental designs 
to be used for the estimation of thermal properties of composite materials. This includes 
both one-dimensional (isotropic) and two-dimensional (anisotropic) analyses. Experiments 
were then conducted for the one-dimensional case, using the optimal design, to estimate 
the effective thermal conductivity perpendicular to the fibers and the effective volumetric 
heat capacity of a composite consisting of IM7 graphite fibers and a Bismaleimide epoxy 
matrix. The estimation procedure used was the modified Box-Kanemasu method. The 
following conclusions can be made based on the obtained results. 

6.1 Optimal Experimental Designs 

In this investigation, optimal experimental designs were determined for both one- 
dimensional and two-dimensional heat conduction processes. In the two-dimensional 
analysis, two different configurations were investigated, both allowing for the estimation 
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of the effective thermal conductivity in two directions and the effective volumetric heat 
capacity. 

6.1.1 One-Dimensional Optimal Experimental Design 

For the one-dimensional experimental design, three experimental parameters were 
optimized: dimensionless heating time, temperature sensor location, and dimensionless 
experimental time. The following conclusions can be made based on the results obtained 
for the specific geometry and boundary conditions used in this analysis: 

1) The optimal dimensionless heating time is 2.2. However, the maximum 
determinant curve had a rather flat peak between heating times of 2.0 and 2.5. 
Therefore, any values within this range can be used. 

2) The optimal temperature sensor location is at the heated surface. 

3) The optimal dimensionless experimental time is approximately 5.0. Note however, 
that this is a conservative choice. 

6.1.2 Two-Dimensional Optimal Experimental Designs 

For the two-dimensional experimental design, two configurations were analyzed. 
Both configurations had a heat flux applied over a portion of one boundary, with the 
remainder of the boundary insulated. In addition. Configuration 1 had known, constant 
temperatures at the remaining three boundaries, while the second configuration had a 
known constant temperature at the boundary opposite to the heat flux and insulated 
conditions at the remaining two boundaries. For each configuration, the optimal 
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experimental parameters determined include the temperature sensor location (x/,y/), the 
dimensionless heating time, the location of the heat flux, and the dimensionless 
experimental time. Based on the obtained results for the specific geometry and boundary 
conditions used in this analysis, the following conclusions can be made: 

6. 1.2.1 Conclusions for Configuration 1 

1) The optimal dimensionless heating time is 1.4. 

2) The optimal temperature sensor location occurs at a x s + of 0.0 (or at the heated 
surface) and ay/ of 0.13 (13% of L y from the bottom edge). 

3) The optimal location of the heat flux is across the entire y + boundary (Lp tl + = 1.0). 

4) The optimal dimensionless experimental time is approximately 4.0. Note however, 
that this is a conservative choice. 

6. 1.2.2 Conclusions for Configuration 2 

1) The optimal dimensionless heating time is 1.55. 

2) The optimal temperature sensor location occurs at a x* of 0.0 (or at the heated 
surface) and a y/ of 0.77 (77% of Ly from the bottom edge). 

3) The optimal location of the heat flux is across 89% of the y + boundary 

(V=0.89). 

4) The optimal dimensionless experimental time is approximately 4.0. Again, this 
is a conservative value. 

The following conclusion can also be made when comparing the two configurations: 
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1) Configuration 1 should be utilized when estimating k x . eff> k y . eff) and C eff 
simultaneously, increasing the accuracy of the resulting property estimates. 

6.2 Thermal Property Estimates 

The estimation of the thermal properties, namely the effective thermal conductivity 
perpendicular to the fibers and the effective volumetric heat capacity, was conducted for 
the one-dimensional analysis using the modified Box-Kanemasu method. This estimation 
procedure requires both measured and calculated temperatures. The measured 
temperatures were obtained from experiments conducted using the optimal experimental 
parameters. The following conclusions can be made based on the results of this portion 
of the investigation: 

1) The effective thermal conductivity perpendicular to fibers (k x .^) is 0.52 W/m°C. 

2) The effective volumetric heat capacity (C eff ) is 1.48 MJ/m 3o C. 

3) The estimated parameters are both reliable, as shown by the small confidence 
intervals and the Root Mean Square values. 

4) The estimates for k x . eff are more accurate than for C eff . 

5) The sequential estimates converge to a steady value, indicating that the heat 

conduction model and experimental design are satisfactory. 

6) No bias error in the calculated temperatures is apparent 
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Chapter 7 


Recommendations 

The estimation procedure used in this investigation to determine the effective 
thermal conductivity perpendicular to the fiber axis and the effective volumetric heat 
capacity proved to be quite accurate. However, because the problem had been simplified 
by conducting the required experiments at room temperature, the actual environmental 
conditions that many composites, especially in aerospace vehicles, are subjected to have 
not been accurately described. These operating conditions usually occur over extreme 
temperature ranges, resulting in temperature dependent thermal properties. Therefore, in 
order to accurately determine the temperature distributions within these structures during 
actual operational conditions, it is necessary to characterize this temperature dependence. 
The estimation procedure can be modified to include this dependence by assuming a 
functional relationship between the thermal properties and temperature. For example, the 
thermal conductivity can be approximated by a piece-wise linear function with 
temperature: 

Kff = k i + (*M - k ) 
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where k { and k i+] are the coefficients to be estimated. The mathematical model can now 
be modified to include temperature dependent properties and the least squares function 
given in Eq. (3.33) can be minimized with respect to k \ and k i+1 . 

In addition, it is also recommended that the results of the two-dimensional 
optimization procedure be verified. This can be done by conducting experiments for the 
two configurations using both the optimal experimental parameters and arbitrary 
experimental parameters. The measured temperatures will then be utilized in the 
estimation procedure to determine the thermal properties, k x .^ k y . eff> and C eff . The 
estimates obtained using the experiments conducted with the optimal parameters should 
provide the smallest confidence intervals. (Recall that the optimization procedure selected 
for this study has the effect of minimizing the confidence intervals of the estimated 
parameters). 

The Box-Kanemasu method could also be implemented to determine the thermal 
properties using both configurations simultaneously. Based on the magnitude of the 
sensitivity coefficients, temperature measurements from Configuration 2 would be used 
to estimate k x . eff and C eff> while k y . eff would be estimated using measurements taken from 
Configuration 1. Using both configurations together will then supply the most accurate 
estimates for all three thermal properties. 

Other properties could also be estimated, such as the thermal contact resistance 
between composite components. Here, the least squares function would be minimized 
with respect to the contact resistance. Furthermore, efforts could be taken for the analysis 
of complex structures. 
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Appendix A 


The FORTRAN program 1DOPT.FOR 


This program, 1DOPT.FOR, is used to calculate the maximum determinant value 


for the one-dimensional analysis, and to determine the corresponding optimal experimental 


parameters. 

PROGRAM ONEDOPT 
C Written by Debbie Moncman, 1993 

DOUBLE PRECISION BETAM,BETA2M, SUMT,SUMC,SUMK 
DOUBLE PRECISION FF1,FF2,BETA2(0: 1000),M,TIME,TTMEH 
DOUBLE PRECISION BETA(0:10(X)),TT,X1T,X2T,INCRX2,PI 
DOUBLE PRECISION Xi,Tl , XI, X2, TTIME, DELTA, INCRET 
DOUBLE PRECISION TMAX,D,XTX1 1,XTX12,XTX22,INCRX1 
COMMON Xi, TIME, TIMER 

OPEN(UNTT = 15, FILE =’T.DAT’, STATUS=’ UNKNOWN’) 

OPEN (UNIT = 20, FILE =’XEDAT\ STATUS**’ UNKNOWN’ ) 
OPEN (UNIT = 25, FILE =’X2.DAT’, STATUS=’UNKNOWN’) 
OPEN(UNIT = 30, FILE =’Dxl.DAT’, STATUS =’ UNKNOWN’) 
OPEN(UNIT = 35, FILE - ’DM AX. DAT’, STATUS=’UNKNOWN’) 
PI =DACOS(-1DO) 

DELTA = 0.0250D0 
TTIME = 6.d0 
Xi = 0.0D+0 

DO 7 TIMEH = DELTA, TTIME, DELTA 
TMAX = 0.0D0 
DMAX = O.ODO 
XTX11 = O.ODO 
XTX12 = O.ODO 
XTX22 = O.ODO 
TT = 0.0D+0 
X1T = O.OD+O 
X2T = 0.0D+0 
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DO 20, TIME = DELTA, TTIME , DELTA 
CALL MODELOT) 

CALL SENS(X1T,X2T) 

XTX11=XTX11+ X1T * X1T 
XTX12=XTX12+ X1T * X2T 
XTX22=XTX22+ X2T * X2T 
C FIND TMAX 

IF(TT.GT.TMAX)TMAX = TT 
D=(l BD0/(TMAX*TMAX*TIME/DELTA))**2* 
+ (XTX 1 1 *XTX22-XIX 12*XTX 12) 

C FIND DMAX 

IF(D.GE.DMAX) THEN 
DMAX = D 
ENDIF 

WRITE(15,12)TIME/DELTA,TIME,TT 
WRJTE(20, 1 2)TIME,X1 T 
WRITE(25,12)TIME,X2T 
WRITE(30,40) TIME,D 
40 FORMAT(lX, D10.5, 4X, D12.5) 

12 FORMAT(2(lX,D12.4),4XJ)12.5) 

14 FORMAT(lXJ)12.4,3(4X,D12.5)) 

20 CONTINUE 

WRTIE(35,8)DMAX, TIMER 
8 FORMAT(lX,2D15.6) 

7 CONTINUE 

CLOSE(15) 

CLOSE(20) 

CLOSE(25) 

CLOSE(30) 

CLOSE(35) 

STOP 

END 


* * * ** * * * * * * * * ♦ * + * * 3|t $ * * * * 9|C * ** * * ** * * * * * £ * * 3|C * * * 3|C * * $ * * * * 4c * * * * * * * * * ))( * * * * * $ 

C Subroutine to calculate the dimensionless temperature 
SUBROUTINE MODEL(TT) 

DOUBLE PRECISION FF1 ,BETA2M,TIME,TIMEH,T1 ,PI 
DOUBLE PRECISION FF2,INCRET,SUMT ,Xi,BETAM,M,TT 
DOUBLE PRECISION BETA2(0:1000), BETA(0:1000) 

COMMON Xi,TIME,TIMEH 
PI =DACOS(-1DO) 

SUMT = 0.0D0 
DO 10, M = 1, 1000, 1 

BETA(M) = (M - 0.5D0)*PI 
BETA2(M) = BETA(M)**2.0D+0 
BET AM = BETA(M) 

BETA2M = BETA2(M) 

FF1 = EXP(-BETA2M * TIME) 

IF (TIME.LE.TIMEH) THEN 
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T1 = FF1 
ELSE 

FF2 = EXP(-BETA2M *(TIME - TIMEH)) 

T1 = FF1-FF2 
ENDIF 

INCRET = Tl*COS(BETAM*Xi)/BETA2M 
IF(M.NE.l) THEN 

IF(ABS(INCRET/SUMT).LT.1.0D-20) then 
GOTO 15 
ENDIF 
ENDIF 

SUMT = SUMT + INCRET 
10 CONTINUE 

15 IF(TIMEJLE.TIMEH) THEN 

TT = 1.0D+0 - Xi - 2.0D+0*SUMT 
ELSE 

TT = -2.0D+0*SUMT 
ENDIF 
RETURN 
END 

£) *** * * * * * * * 4c * * * * ** *** * * **** ** ** * * *********** * ********* ** ****** ** 3|C 9|C s|c 3(c>(e >|c 

C Subroutine to Calculate the dimensionless Sensitivity Coefficients 
SUBROUTINE SENS(X1T,X2T) 

DOUBLE PRECISION TIME,TIMEH,BETAM,BETA2M, SUMK,Xi,M 
DOUBLE PRECISION SUMC,FF1,FF2,X1,X2,INCRX1,INCRX2 
DOUBLE PRECISION BETA(0:1000),BETA2(0:1000),PI 
DOUBLE PRECISION X1T, X2T 
COMMON Xi, TIME, TIMEH 
PI =D ACOS (- 1 DO) 

SUMK = 0.0D+0 
SUMC = 0.0D+0 
DO 30, M = 1, 1000, 1 

BETA(M) = (M - 0.5D0)*PI 
BETA2(M) = BETA(M)**2.0D+0 
BETAM = BETA(M) 

BETA2M = BETA2(M) 

FF1 = EXP(-BETA2M * TIME) 

IF(TIME.LE.TIMEH) THEN 
XI = FF1 * (1/BET A2M + TIME) 

X2 = TTME*FF1 
ELSE 

FF2 = EXP(-BETA2M *(TIME - TIMEH)) 

XI = ( l/BETA2M+TIME)*I 7 Fi -( 1/B ET A2M+(TIME-TIMEH))*EE2 

X2 = TIME*FF1 - (TIME-TIMEH)*FF2 

ENDIF 

INCRX1 = XI *COS(BETAM * Xi) 

INCRX2 = X2*COS(BETAM * Xi) 

IF(M.NE.l) THEN 
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IF(ABS(INCRX1/SUMK).LT.1.0D-20.AND. 
+ AB S (INCRX2/SUMC).LT. 1 .0D-20) THEN 

GOTO 16 
END IF 
END IF 

SUMK = SUMK + INCRX1 
SUMC = SUMC + INCRX2 
30 CONTINUE 
16 IFCnMElE.lTMEH) THEN 

X1T = -(1.0D0 - Xi) + 2.0D0*SUMK 
X2T = -2.0D0 * SUMC 
ELSE 

X1T = 2.0D0*SUMK 
X2T = -2.0D0 * SUMC 
END IF 
RETURN 
END 
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Appendix B 


The FORTRAN program 2DC10PT.F0R 


This program, 2DC10PT.FOR, is used to calculate the maximum determinant 
value for Configuration 1 of the two-dimensional analysis, and to determine the 
corresponding optimal experimental parameters. 


PROGRAM CFONEOPT 
C Written by Debbie Moncman, 1994 

DOUBLE PRECISION PI,SUMC,SUMKX,SUMKY 

DOUBLE PRECISION Lr,K,BETAN,BETAN2JEXPON,TERMl,TX 

DOUBLE PRECISION FFl^KX,XKY^CC,FF2,CONST,INCRKX,INCRKY,INCRC 

DOUBLE PRECISION TERM2,TERM3 ,S UMT,TEMP,INCRT 

DOUBLE PRECISION EXPONTM,EXPONTH,SSSUMT, SSINCRT 

DOUBLE PRECISION SSSUMKX,SSINCRKX,SSSUMKY,SSINCRKY 

DOUBLE PRECISION TERM4 t X,Y,Yp f TERM5,TMAX 

DOUBLE PRECISION X1T^2T^C3T,XTX11,XTX12^TX13^TX22 

DOUBLE PRECISION XTX23,XTX33,DETADMAX,THOPT 

DOUBLE PRECISION TIME/IIMEH/ITME.DELTA, TIMET 

INTEGER M,N 

C OPEN(UNIT==40 JTLE=’ TCTL1K17 DAT ,STATUS=’UNKNOWN’ ) 

C OPEN(UNrr=654TLE=’dDAT’,STATUS=’UNKNOWN’) 

PI = DACOS(-l.ODO) 

DELTA = 0.05D0 
TOME = 6.0D0 
K = 7.0d0 
Lr = 0.048276170 
SSSUMT = O.IX) 

DO 2, M = 1, 3000 
Y = 0.5170 
Yp = 1.0D0 
X = 0J70 

TERM1 = DSIN(M*PI*Y) 

TERM2 = EDO - DCOS(M*PI*Yp) 
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IF(TERM1 .EQ.O..OR.TERM2.EQ.O.) GOTO 2 
TERM3 = M**2*PI**2*Lr**2*K 
DO 3, N = 1, 3000 

BET AN = PI*(N-0.5D0) 

BETAN2 = BETAN*BETAN 
TERM4 = TERM3 + BETAN2 
TERM 5 = DCOS(BETAN*X) 

IF(TERM5JEQ.O.) GOTO 3 

S SINCRT=TERM1 *TERM2*TERM5 * ( 1 .D0/(M*TERM4)) 
SSSUMT = SSSUMT + SSINCRT 
3 CONTINUE 

2 CONTINUE 

TMAX = SSSUMT*(4.D0/PI) 

DO 150, Yp = 0.1D0, 1.0D0, 0.05D0 
DO 125, X = 0.0D0, 1.0D0, 0.05DO 
DO 100, Y = 0.0D0, 1.0D0, 0.05D0 
DMAX = OJX) 

SSSUMT = 0.0D0 
SSSUMKX = 0.0D0 
SSSUMKY = 0.D0 
DO 400, M = 1, 3000 
TERM1 = DSIN(M*PI*Y) 

TERM2 = EDO - DCOS(M*PI*Yp) 

IF(TERM1JEQ,O.OR.TERM2EQ.O.) GOTO 400 
TERM3 = M**2*PI**2*Lr**2*K 
DO 300, N = 1, 3000 
BETAN = PI*(N-0.5D0) 

BETAN2 = BETAN*BETAN 
TERM4 = TERM3 + BETAN2 
TERM5 = DCOS (BETAN*X) 

IF(TERM5.EQ.O.) GOTO 300 
SSINCRT=TERM1 *TERM2*TERM5*(1 D0/(M*TERM4)) 

SSINCRKX = TERM l*TERM2*DCOS(BETAN*X)*( 1 X)0/(M*TERM4)) 
+ *((TERM3/TERM4)-1 .DO) 

SSINCRKY = TERM 1 * TERM2 *DCOS (BETAN*X) *( 1 JD0/(M*TERM4)) 

+ *(-TERM3/TERM4) 

SSSUMT = SSSUMT + SSINCRT 

SSSUMKX = SSSUMKX + SSINCRKX 
SSSUMKY = SSSUMKY + SSINCRKY 
300 CONTINUE 
400 CONTINUE 

DO 650, TIMER = DELTA, TOME, DELTA 
XTX11 = OOO 
XTX12 = OJX) 

XTX13 = OJX) 

XTX22 = OJX) 

XTX23 = OJX) 

XTX33 = OJX) 

SUMT = 0.0D0 
SUMC = OJX) 
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SIJMKX = O.DO 
SUMKY = O.DO 

DO 200, TIME = DELTA, TTIME, DELTA 
DO 500, M = 1, 10000,1 

TERM! = DSIN(M*PI*Y) 

TERM2 = 1.D0 - DCOS(M*PI*Yp) 

IF(TERM1.EQ.O..OR.TERM2.EQ.O.)GOTO 500 
TERM3 = M**2*PI**2*Lr**2*K 
DO 600, N = 1, 1000, 1 

BETAN = PI*(N-0.5D0) 

BETAN2 = BETAN+BETAN 
TERM4 = TERMS + BETAN2 
TERM5 = DCOS(BETAN*X) 

IF(TERM5.EQ.O.)GOTO 600 
EXPON = TERM4 
EXPONTM = EXP0N*T1ME 
IF(TIMEI J E.TIMEH) THEN 

IF(EXPONTMJLT.225.) THEN 
FF1 = DEXP(-EXPONTM) 

ELSE 
FF1 = OJDO 
END IF 
TX = FF1 

XKX = BETAN2*TIME*FFl-((TERM3/EXPON)-lJX))*FFl 
XKY = TERM3 *TIME*FF1 + ((TERM3/EXPON)*FFl) 

XC = -EXPON*TIME*FFl 

ELSE 

EXPONTH = (EXPON* (TIME-TIMEH)) 

IF(EXPONTHET.225..AND£XPONTMiT.225.)THEN 
FF1 = DEXP(-EXPONTM) 

FF2 = DEXP(-EXPONTII) 

ELSE IF(EXPONTH.LT.225. . AND.EXPONTM.GE.225 .) THEN 
FF2 = DEXP(-EXPONTH) 

FF1 = OJDO 

ELSE IF(EXPONTH.GE.225..AND.EXPONTM.LT.225.) THEN 
FF1 = DEXP(-EXPONTM) 

FF2 = OJDO 
ELSE 
FF1 = 0.D0 
FF2 = O.DO 
ENDIF 

TX = FF2 - FF1 

XKX = ((TERM3/EXPON)-1.DO)*(FF2-FF1) + 

+ BETAN2*TIME*FF1 - BETAN2*(TIME-TIMEH)*FF2 

XKY = (-TERM3/EXP0N)*(FF2-FF1) + TERM3 *TIME*FF1 
+ - FERM3 * (TIME-TIMEI I)*FF2 

XC = EXPON* (TME-HMEI I) *FF2 - EXPON*ITME*FFl 
ENDIF 

CONST = TERM1 *TERM2*TERM5 *(1 ,D0/(M*TERM4)) 

INCRT = TX*CONST 
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INCRC = XC*CONST 
INCRKX = XKX* CONST 
INCRKY = XKY*CONST 

IF(SUMT.NE.O. .AND. SUMKX.NE.0. . AND. SUMKY.NE.O. .AND . 

+ SUMC.NE.O.)THEN 

IF(ABS(INCRT/SUMD.LT.1X)-20.AND.ABS(INCRKX/SUMKX)XT. 
+ 1X)-20.AND.ABS(INCRKY/SUMKY)JLT.1.D-20.AND.ABS 

+ (INCRC/SUMC)lT.lX)-20) THEN 

GO TO 410 
ENDIF 
ENDIF 

SUMT = SUMT + INCRT 
SUMC = SUMC + INCRC 
SUMKX = SUMKX + INCRKX 
SUMKY = SUMKY + INCRKY 
600 CONTINUE 
410 IF(N.EQ.1)THEN 

IF(AB S (INCRKX) XT. 1 .D-20. AND . AB S (INCRKY) ET. 1 .D-20 
+ . AND. AB S (INCRC).LT. 1 J)-20 . AND. AB S (INCRT) ITT. 

+ 1X)-20)THEN 

GO TO 450 
ENDIF 
ENDIF 

500 CONTINUE 

450 IF(TIME.LE.TIMEH) THEN 

TEMP = (4.D0/PI)*(SSSUMT-SUMT) 

X3T = (4.0D0/PI)*SUMC 

X1T = (4XH)/PI)*(SSSUMKX + SUMKX) 

X2T = (4JX)/PI)*(SUMKY + SSSUMKY) 

ELSE 

TEMP = (4 .0D0/PI) * S UMT 
X3T = (4D0/PI)*SUMC 
X1T = (4D0/PI)*SUMKX 
X2T = (4.D0/PI)*SUMKY 
ENDIF 

WRITE(40, 14)y,TIME,TEMP,xlt,x2t,x3t 
14 FORMAT(lx,f5.2,5el3.5) 

X1T = X1T/TMAX 
X2T = X2T/TMAX 
X3T = X3T/TMAX 
XTX11 = XTX11 + X1T*X1T 
XTX12 = XTX12 + X1T*X2T 
XTX13 = XTX13 + X1T*X3T 
XTX22 = XTX22 + X2T*X2T 
XTX23 = XTX23 + X2T*X3T 
XTX33 = XTX33 + X3T*X3T 

DET = XTX1 1*(XTX22*XTX33 - XTX23*XTX23) - XTX12*(XTX12*XTX33 
+ - XTX13*XTX23) + XIX13*(XTX12*XTX23 - XTX13*XTX22) 

D = (1 .D0/(TIMF7DELTA))**3*DET 
IF(D.GE.DMAX) THEN 
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DMAX = D 
THOPT = TIMER 
TIMET = TIME 
ENDIF 

SUMT = 0.0D0 
SUMC = 0.D0 
SUMKX = 0.D0 
SUMKY = 0.D0 
200 CONTINUE 
650 CONTINUE 

WRITE(65,1 10)X,Y,Yp,DMAX,THOPT,TIMET 

110 FORMAT(lX,3(2X,F6.3),3E13.6) 

100 CONTINUE 
125 CONTINUE 
150 CONTINUE 
STOP 
END 



Appendix C 


The FORTRAN program 2DC20PT.F0R 


This program, 2DC20PT.F0R, is used to calculate the maximum determinant 
value for Configuration 2 of the two-dimensional analysis, and to determine the 


corresponding optimal experimental parameters. 

PROGRAM CFTWOOPT 
C Written by Debbie Moncman, 1994 

DOUBLE PRECISION PI,SUMC,SUMKX,SUMKY 

DOUBLE PRECISION Lr,K,BETAN,BETAN2,TERMl,Tl 

DOUBLE PRECISION FF1 r XKX,XKY ,XCP, FF2, CONST, INCRKX,INCRKY,INCRC 

DOUBLE PRECISION TERM2,TERM3,SUMT,TEMP,INCRT,KYN,CPN 

DOUBLE PRECISION EXPONTM,EXPONTI I,S S S IJMT, SSINCRT 

DOUBLE PRECISION SSSUMKX,SSINCRKX,SSSUMKY,SSINCRKY,SSTEMP2 

DOUBLE PRECISION TERM4,X,Y,Xp,TMAX,TEMPl,TEMP2,SSTEMPl 

DOUBLE PRECISION X1T,X2T,X3T,XTX1 1,XTX12,XTX1 3,XTX22 

DOUBLE PRECISION XTX23,XTX33,DET,D,DMAX,THOPT,SSUMTN 

DOUBLE PRECISION TIME, TIMEH/IIIME, DELTA, timet, term5 

INTEGER M,N 

OPEN(UNIT=40JFILEs=’TILlK17J)AT’,STATUS=’UNKNOWN’) 
OPEN(UNrr=65,Fn.E=’DIVARYY.DAT’,STATUS= ’UNKNOWN’) 

PI = DACOS(-l.ODO) 

DELTA = 0.05D0 
ITIME = 6.0DO 
K = 7.0d0 
Lr = 0.048276D0 
SSSUMT = 0.D0 
SSUMTN = 0.D0 
DO N = 1, 3000,1 
Y = 0.0D0 
Xp = 1.0D0 
X = 0.5D0 
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BETAN = PI*(N-0.5D0) 

BETAN2 = BETAN*BETAN 
INCRT = DC0S(BETAN*Y)/BETAN2 
SSUMTN = SSUMTN + INCRT 
ENDDO 

SSTEMP1 = 2JD0*Xp*SSUMTN 
DO 2, M = 1, 3000 
Y = 0.0D0 
Xp = 1.0D0 
X = 0.5D0 

TERM1 = DCOS(M*PI*X) 

TERM2 = DSIN(M*PI*Xp) 

IF(TERM 1 .EQ.O..OR.TERM2.EQ.O.) GOTO 2 
TERM3 = M**2*PI**2*Lr**2*K 
DO 3, N = 1, 3000, 1 

BETAN = PI*(N-0.5D0) 

BETAN2 = BETAN*BETAN 
TERM4 = TERM3 + BETAN2 
TERM5 = DCOS(BETAN*Y) 

IF(TERM5.EQ.O.) GOTO 3 

SS INCRT = TERM 1 *TERM2 * term5 *( 1 ,D0/(M*TERM4)) 
SSSUMT - SSSUMT + SSINCRT 
3 CONTINUE 

2 CONTINUE 

SSTEMP2 = 4.0D0/PI*SSSUMT 
TMAX = SSTEMP2 + SSTEMP1 
DO 150, Xp = 0.1D0, 1.05D0, 0.05D0 
DO 125, Y = 0.0D0, 1.0D0, 0.05DO 
DO 100, X = 0.0D0, 1.0D0, 0.05D0 
DMAX = 0J30 
SSSUMT = 0.0D0 
SSSUMKX = O.ODO 
SSSUMKY = 0.D0 
IX) 400, M = 1, 3000 
TERM1 = DCOS(M*PI*X) 

TERM2 = DSIN(M*PI*Xp) 

IF(TERMEEQ.O..OR.TERM2.EQ.O.) GOTO 400 
TERM3 = M**2*PI**2*Lr**2*K 
DO 300, N = 1, 3000, 1 
BETAN = PI*(N-0.5DO) 

BETAN2 = BETAN* BETAN 
TERM4 = TERM3 + BETAN2 
TERM 5 = DCOS (BETAN* Y) 

IF(TERM5 JBQ.O.) GOTO 300 

SSINCRT = TERM! *TERM2*term5*(l .D0/(M*TERM4» 
SS1NCRKX = TERM 1 *TERM2*term5 *( 1 .D0/(M*TERM4» 
+ *(-TERM3/TERM4) 

SSINCRKY = TERM 1 *TERM2*teim5 *( 1 .D0/(M*TERM4» 
+ *((TERM3/TERM4)-1.D0) 

SSSUMKX = SSSUMKX + SSINCRKX 
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SSSUMKY = SSSUMKY + SSINCRKY 
SSSUMT = SSSUMT + SSINCRT 

300 CONTINUE 
400 CONTINUE 

DO 650, TIME! I = DELTA, TOME, DELTA 

XTX11 = 0.D0 

XTX12 = 0D0 

XTX13 = 0D0 

XTX22 = 0 DO 

XTX23 = 0D0 

XTX33 = 0D0 

SUMT = O.ODO 

SUMC = O.DO 

SUMKX = O.DO 

SUMKY = ODO 

DO 200, TIME = DELTA, TOME, DELTA 
DO N = 1, 500,1 
BETAN = PI*(N-0.5D0) 

BETAN2 = BETAN* BETAN 
EXPONTM = BETAN2*TIME 
IF(TIME.LE.TIMEH) THEN 

rF(EXPONTM.LE.225.) THEN 
FF1 = DEXP(-EXPONTM) 

ELSE 

FF1 = O.DO 
ENDIF 

T1 = (1 DO - FF1) 

XCP = TTME*FF1 

XKY = (-1D0/BETAN2) + (TIME +(1.D0/BETAN2))*FF1 
ELSE 

EXPONTH = B ETAN2 * (TIME-TIMEH) 

IF(EXPONTHTE.225..AND£XPONTMJLE.225.) THEN 
FF1 = DEXP(-EXPONTM) 

FF2 = DEXP(-EXPONTH) 

ELSE I FCEXTONTI I . GT. 22 5 . . AND .EXPONTM .LE. 225 . )TI I IiN 
FF1 = DEXP(-EXPONTM) 

FF2 = ODO 

ELSE IF(EXPONTH.LE.225..AND.EXPONTM.GT.225.) THEN 
FF1 = O.DO 

FF2 = DEXP(-EXPONTH) 

ELSE IF(EXTONTH.GT.225..AND.EXPONTM.GT.225.) THEN 
FF1 = ODO 
FF2 = ODO 
ENDIF 

T1 = FF2 - FF1 

XKY = ((1D0/BETAN2)+TIME)*FF1 - ((DME-TIMEH) 

+ + (1D0/BETAN2))*FF2 

XCP = TIME*FF1 - CUME-TIMEI I) *FF2 
ENDIF 

INCRT = T1 *DCOS(BETAN* Y)/BETAN2 
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INCRKY = DCOS (BETAN* Y)*XKY 
INCRC = XCP*DCOS(BETAN*Y) 

IF(SUMTJSfE.O..AND.SUMKY.NE.O..AND.SUMCJ^E.O.)THEN 

IF(ABS(INCRT/SUMT).LE.1.D-10.AND.ABS(INCRKY/SUMKY) 
.LE.1X)-10.AND.ABS(INCRC/SUMC)X,E.1.D-10) THEN 
GOTO 13 
ENDIF 
ENDIF 

SUMT = SUMT + INCRT 
SUMKY = SUMKY + INCRKY 
SIJMC = SUMC + INCRC 
ENDDO 

TEMPI = 2X)0*Xp*SUMT 
KYN = 2D0*Xp*SUMKY 
CPN = -2.D0*Xp*SUMC 
SUMT = 0.D0 

sumc = o do 

SUMKY = 0.D0 
DO 500, M = 1, 10000,1 
TERM1 = DCOS(M*PI*X) 

TERM2 = DSIN(M*PI*Xp) 

IF(TERM1£Q.O..OR.TERM2JEQ.O.) GOTO 500 
TERM3 = M**2*PI**2*Lr**2*K 
DO 600, N = 1, 1000, 1 

BETAN = PI*(N-0.5D0) 

BETAN2 = BETAN*BETAN 
TERM4 = TERM3 + BETAN2 
TERM5 = DCOS (BET AN* Y) 

IF(TERM5£Q.0.) GOTO 600 
EXPONTM = TERM4*TIME 
IFCTIME.LE.TIMEH) THEN 

IF(EXPONTMJLT.225.) THEN 
FF1 = DEXP(-EXPONTM) 

ELSE 
FF1 = 0JX) 

ENDIF 
T1 = FF1 

XKX = (TERM3/TERM4)*FF 1 + TERM3 *TIME*FF1 
XKY = ((TERM3/TERM4)-1D0)*(-FF1)+BETAN2*TIME*FF1 
XCP = -(TERM4*TIME*FF1) 

ELSE 

EXPONTH = TERM4*(TIME-TIMEH) 

Il 7 (EXPONTH.LE .225 . . AND .EXPONTM .LE.225 .) THEN 
FF1 = DEXP(-EXPONTM) 

FF2 = DEXP(-EXPONTH) 

ELSE IF(EXPONTH.GT.225 .. ANDEXPONTMiE.225 .)THEN 
FF2 = 0.D0 

FF1 = DEXP(-EXPONTM) 

ELSE IF(EXPONTH.LE.225..AND.EXPONTM.GT.225.) THEN 
FF1 = 0.D0 



FF2 = DEXP(-EXPONTH) 

ELSE IF(EXPONTH.GT225..ANDHXPONTM.GT.225.) THEN 
FF1 = OXK) 

FF2 = 0XH) 

ENDIF 

T1 = FF2 - FF1 

XKX = (-TERM3/TERM4)*(FF2-FF1) + TERM3 *TIME*FE1 
+ - TERM3*(TIME-TIMEH)*FF2 

XKY = ((TERM3/TERM4)-1.D0)*(FF2-FF1) + BETAN2*TIME*FF1 
+ - BETAN2*(TTME-TIMEH)*FF2 

XCP = TERM4*(TIME-TIMEH)*FF2 - TERM4*TTME*FF1 
ENDIF 

CONST = TERM1 *TERM2*TERM5*(1 .D0/(M*TERM4» 

INCRT = Tl*CONST 
INCRKX = XKX*CONST 
INCRKY = XKY*CONST 
INCRC = XCP*CONST 

IF(SUMT.NEB..AND.SUMKXT4E.0..AND.SUMKY.NE.0..AND. 

+ SUMC.NE.O.)THEN 

IF(AB S(INCRT/SUMT)XT. 1 .D-20.AND.ABS(INCRKX/SUMKX)XT. 
+ l.D-20.AND.ABS(INCRKY/SUMKY)XT.l.D-20.AND.ABS 

+ (INCRC/SUMC)XT. 1 X)-20) THEN 

GO TO 410 
ENDIF 
ENDIF 

SUMT = SUMT + INCRT 
SUMC = SUMC + INCRC 
SUMKX = SUMKX + INCRKX 
SUMKY = SUMKY + INCRKY 
600 CONTINUE 
410 IF(N.EQ.1)THEN 

IF(AB S (INCRKX)XT. 1 .D-20. AND. AB S(INCRKY)XT. 1 .D-20 
+ .AND.ABS(INCRC)XT.1X)-20.AND.ABS(INCRT)XT. 

+ 1X)-20)THEN 

GO TO 450 
ENDIF 
ENDIF 

500 CONTINUE 

450 IF(TIMEXE.TIMEH) THEN 

TEMP2 = (4.D0/PI)*(SSSUMT - SUMT) 

X1T = (4X0/PI)*(SSSUMKX + SUMKX) 

X2T = KYN + (4J)0/PI)*(SS SUMKY + SUMKY) 

X3T = CPN + (4.D0/PI)*SUMC 
ELSE 

TEMP2 = (4.D0/PI)*SUMT 
X1T = (4X)0/PI)*SUMKX 
X2T = KYN + (4X)0/PI)*SUMKY 
X3T = CPN + (4.D0/PI)*SUMC 
ENDIF 

TEMP = TEMPI + TEMP2 
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WRITE(40, 14)y,TIME,TEMP,x It,x2t,x3t 
14 format(lx,f5.2,5el3.5) 

X1T = X1T/TMAX 
X2T = X2T/TMAX 
X3T = X3T/TMAX 
XTX11 = XTX11 + X1T*X1T 
XTX12 = XTX12 + X1T*X2T 
XTX13 = XTX13 + X1T*X3T 
XTX22 = XTX22 + X2T*X2T 
XTX23 = XTX23 + X2T*X3T 
XTX33 = XTX33 + X3T*X3T 

DET = XTX11*(XTX22*XTX33 - XTX23*XIX23) - XTX12*(XTX12*XTX33 
+ - XTX13*XTX23) + XTX13*(XTX12*XTX23 - XTX13*XTX22) 

D = (1JX)/(TIME/DELTA))**3*DET 
IF(D.GE.DMAX) THEN 
DM AX = D 
THOPT = TTMEH 
TIMET = TIME 
ENDIF 

SUMT = 0.0D0 
SIJMC = 0.IX) 

SUMKX = 0.D0 
SUMKY = 0.D0 
200 CONTINUE 
650 CONTINUE 

WRITE(65,1 10)Xp,Y^X,DMAX,THOPT,TIMET 
110 FORMAT(lX,3(2X r F6.3),3E13.6) 

100 CONTINUE 
125 CONTINUE 
150 CONTINUE 
STOP 
END 
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Appendix D 


The FORTRAN program MODBOX.FOR 


This program, MODBOX.FOR, uses the modified Box-Kanemasu method to 
estimate the thermal properties. 


PROGRAM NLINA 

CCCCCCCCC PROGRAM DESCRIPTION CCCCCCCCC 

C THIS PROGRAM USES THE MATRIX INVERSION LEMMA (BASED ON C 
C THE GAUSS LINEARIZATION METHOD) AND THE BOX-KANEMASU C 
C METHOD TO ESTIMATE THE PARAMETERS OF A GIVEN MODEL. C 


£****************************************************************£ 


C 

C Written by Debbie Moncman, 1993 
C Based on the program, NLINAFOR, by J. V. Beck (1993) 

CCCCCCCCC DIMENSION BLOCK CCCCCCCCC 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION T(1500,5),Y(1500),SIG2(1500),B(5),SC(5),A(5),BS(5), 

+ VINV(5,5),BSS(5),SUMG(5),R(5,5) r EXTRA(20),ERR(1500), 

+ PS(5,5),P(5,5),XTX(5,5),XTY(5),SUM(5),VALUEK(5),BSV(5) 

C II ARACTER*40 INFILE, OI JTFIL 




C COMMON BLOCK C 

COMMON T^,SC,BS,I,ETA,PS,P,B,A,Y,SIG2MODL,VINV,NPJEXl’RA 

COMMON/ERROR/ERR 

COMMON/MOD/ A A,TL, SUM 




C DATA BLOCK C 

DATA EPSDEN,CRITER/1.0D-30,0.0001D+0/ 
£******************************* *******************************,(<*£ 


C INITIALIZATION BLOCK C 

WRITE(*,*)’ENTER THE NAME OF THE INPUT DATA FILE’ 

READ( * ’ ( A40) ’ )INFILE 

OPEN(8,I 7 ELE=INFILE) 
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WRITE(*,*)’ENTER THE NAME OF THE OUTPUT FILE’ 

READ(V(A40)’)OUTFIL 

OPEN(7 JTLE=OUTFIL) 


£****************************************************************£ 


C PROCESS BLOCK C 

C - START READING INPUT VALUES 
C BLOCK 1 

WRITE(7,*)TNPIJT QUANTITIES’ 

READ(8,*) N, NP, NI, MAXIT, MODL, IPRINT 
WRITE(7,*) 

WRTTE(7,*)’BLOCK 1’ 

WRITE(7,*) 

WRTTE(7,*)’N - NUMBER OF DATA POINTS (MEASUREMENTS)’ 
WRITE(7,*)’NP - NUMBER OF PARAMETERS’ 

WRITE(7,*)’NI - NUMBER OF INDEPENDENT VARIABLES’ 
WRITE(7,*)’MAXIT - MAXIMUM NUMBER OF ITERATIONS’ 


WRITE(7,*)’MODL - MODEL NO.(NEEDED IF SEVERAL MODELS ARE USED)’ 
WRITE(7,*)TPRINT - 1 FOR USUAL PRINTOUTS, 0 FOR LESS’ 

WRITE(7,*) 

IF(N.LE.O)THEN 
STOP 
END IF 


WRITE(7,’(/,9X,”N”,8X,”NP”,8X,”NI”,5X,”MAXIT”,5X, 

+”MODL”,4X, ’’IPRINT”)’) 

WRITE(7,’ (7I10)’)N,NP,NI, MAXIT, MODL, IPRINT 
C BLOCK 2 (INITIAL PARAMETER ESTIMATES) 

WRITE(7,*) 

WRITE(7,*)’ BLOCK 2’ 

WRITE(7,*) 

WRITE(7,*)’B(1),...,B(NP) ARE THE INITIAL PARAMETER ESTIMATES.’ 
WRITE(7,*) 

READ(8,*)(B(I),I=1,NP) 

WRITE(7,’(10X,”B(”,I1,”) =”J 7 16.5)’)(I,B(I),I-1,NP) 

C BLOCK 3 (INPUT MEASUREMENTS) 

WRITER,*) 

WRITE(7,*)’BLOCK 3’ 

WRITE(7,*) 

WRITE(7,*)’J - DATA POINT INDEX’ 

WRITE(7,*)’Y(J) - MEASURED TEMPERATURE VALUE’ 
WRTIE(7,*)’SIGMA(J) - STANDARD DEVIATION OF Y(J)’ 

WR1TE(7, *)’ T(J, 1) - FIRST INDEPENDENT VARIABLE’ 

WRITE(7,*) 

WRITE(7,’(/,9X,”J”,6X,”Y(J)”,3X,”SIGMA(J)”,6X, 

+ ”T(J,1)”,6X,”T(J,2)”)’) 

DO 11=1, N 

READ(8,*)J,Y(J),SIG2(J),(T(J,KT),KT=1,NI) 

WRTIE(7,’(I10,7F10.5)’)J,Y(J),SIG2(J),(T(J,KT),KT=1,NI) 

SIG2(J) = SIG2(J)*SIG2(J) 

END DO 

C BLOCK 4 (INPUT ANY EXTRA DATA NEEDED IN THE MODEL) 
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READ(8,*)IEXTRA 

C I EXTRA is the number of constants used in the model such as 

C initial temperatures, surface temperatures, or a heat flux. 

C It equals 0 for no extra input 
WRITE(7,*) 

WRITE(7,*)’ BLOCK 4’ 

WRITE(7,*) 

WRITE(7,*)’IEXTRA - NUMBER OF EXTRA® PARAMETERS (0 IF NONE)’ 
WRITE(7,*) 

WRITE(7,’(10X,”IEXTRA =”,U0)’)IEXTRA 
IF(IEXTRA.GE. 1) THEN 
WRITE(7,*) 

WRTTE(7,*)’EXTRA(1),... ARE EXTRA CONSTANTS USED IN THE MODEL’ 
WRITE(7,*) 

READ(8, *)(EXTRA(I),I= 1 , IEXTRA) 

WRITE(7,’(”EXTRA(”,I2,”) =”,F16.5)’) 

+ (I,EXTRA(I),I= 1 , IEXTRA) 

ENDIF 


C End input begin calculations 
WRITE(7,*) 

WRTTE(7,*)’END INPUT QUANTITIES, BEGIN OUTPUT CALCULATIONS’ 
WRTTE(7,*) 

WRITE(7,*)’SSY - SUM OF SQUARES FOR PRESENT PARAMETER VALUES’ 
WRITE(7,*)’SSYP - SUM OF SQUARES FOR BOX-KANEMASU PARAMETER VAL.’ 
WRITE(7,*)’ SSYP DECREASES TOWARDS A POSITIVE CONSTANT’ 
WRITE(7,*)’ AND SHOULD BE LESS THAN SSY’ 

WRITE(7,*)’G - MEASURE OF THE SLOPE, IT SHOULD APPROACH ZERO’ 
WRITE(7,*)’ AT CONVERGENCE’ 

WRITE(7,*)’H - SCALAR INTERPOLATION FACTOR; ITS A FRACTION OF’ 
WRITE(7,*)’ THE GAUSS STEP GIVEN BY THE BOX-KANEMASU METHOD’ 
WRITE(7,*) 

WRITE(7,*) 

CCCCCCCC PART I OF PROGRAM (GAUSS METHOD) CCCCCCCCC 

CCCCCCCC**** **************************************** *****ccccccccc 
C - Set the P matrix equal to zero 
DO I2=1,NP 
DO K2-1,NP 
PS(K2,I2)=0 
P(K2,I2)=0 
ENDDO 
ENDDO 


DO I3=1,NP 

PS(I3,I3)=B(I3)*B(I3)*L0D+7 
ENDDO 
DOK= 1, NP 
BS(K)=B(K) 

BSV(K)=BS(K) 

SUMG(K) = 0.013+0 
ENDDO 

C — Set XTX and XTY sums equal to 0 
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DOK= 1, NP 
XTY(K) - O.OD+O 
IX) J = 1, NP 
XTX(J,K) = O.OD+O 
ENDDO 
ENDDO 


C — 1 and MAX are counters 
1 = 0 
MAX = 0 

100 MAX = MAX + 1 
SSY = 0.0D+0 
DO 13 = 1, N 
1=13 

CALL MODEL 
CALL SENS 
RES ID = Y(I) - ETA 
SSY = SSY + RES ID * RES ID/S IG2(I) 

C — Calculate XTX, XTY, and SUMG (used in the Box-Kanemasu method) 
DO K = 1, NP 

XTY(K) = XTY (K) + SC(K)*RESID/SIG2(D 
SUMG(K) = SUMG(K) + SC(K)*RESID/SIG2(I) 

IX) L =1, NP 

XTX(K,L) = XTX(K,L) + SC(K)*SC(L)/SIG2(I) 

ENDDO 
ENDDO 
DO K = 1, NP 
A(K) = 0.0D+0 
ENDDO 

C — Calculate ’A’ used in the MIL method 
1X3 K = 1, NP 
DO J = 1, NP 

A(K) = A(K) + SC(J)*P(KJ) 

ENDDO 

ENDDO 

DELSUM = 0.0D+0 

C — Calculate ’DELTA* used in the MIL method 
DO K = 1JNP 

DELSUM = DELSUM + SC(K)*A(K) 

ENDDO 

DELTA = SIG2(I) + DELSUM 
C — Calculate ’K’ used in the MIL method 
DOK= 1, NP 

VALUEK(K) = A(K)/DELTA 
ENDDO 

SUMH = 0.0D+0 

C - Calculate ’SUMH’ used in ’HU’ 

DO J = 1, NP 

SUMH = SUMH + SC(J)*(B(J) - BS(J)) 

ENDDO 
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HU = Y(I) - ETA - SUMH 

C -- Estimated parameters found using the Gauss Method: 

DO K = 1, NP 

B(K) = B(K) + V ALUEK(K)*HU 
ENDDO 

C — Calculate the new P matrix 
DO U = 1,NP 
DO V= 1,NP 

P(U,V) = PS(U,V) - VALUEK(U)*A(V) 

ENDDO 

ENDDO 

C — Make the P matrix symmetrical 
DO J = 2, NP 
JK - J - 1 
DO K = 1, JK 
P(K,J)=P(J,K) 

ENDDO 
ENDDO 
DO J = 1JSTP 
DO K = 1,NP 
PS(KJ) = P(KJ) 

ENDDO 

ENDDO 

£ * * * * * * * * * * * * ♦ * * * * * * * * * 4c * * £ * * * * 4c * j|c * $ * $ ** * * * sft $ $ * * * * * * * * * * afe * * * * * * 9|e * afc * i|c * £ % * 


C -- Done with Gauss calculations, Print results 
IF(IPRINT.GT.O) THEN 
IF(IJEQ.l) THEN 
WRITE(7,*) 

WRITE(7, *)’ SEQUENTIAL ESTS. OF THE PARAMETERS GIVEN BELOW’ 

W RITE(7, *) ’ (THESE EST. ARE FOUND USING THE GAUSS METHOD)’ 
WRITE(7,’ (//,3X,’ T’,6X,’ ’ETA’ ’ ,5X,” RESIDUALS ’ ’ ,7X, 

+ ”B(1)”,8X,”B(2)”,6X,”SC(1)”, 6X,”SC(2)”)’) 

END IF 

WRITE(7,’(I4,6E12.5)’)I, ETA, RESID, (B(IP),IP=1,NP),SC(1)*B(1), 

+ SC(2)*B(2) 

ENDIF 

ENDDO 

WRITE(7,*) 

WRITE(*,*)’END BASIC LOOP’ 

WRITE(7,*)’THE FINAL SEQUENTIAL ESTIMATES WILL NOW BE MODIFIED 
+ USING THE BOX-KANEMASU METHOD’ 

CCCCCCCCC PART II: BOX-KANEMASU MODIFICATION CCCCCCCCC 


CCCCCCCCC** ******************************************** **CCCCCCCCC 


C — Set BSS equal to the initial estimate for that iteration 
DO J = 1, NP 
BSS(J) = BS(J) 

ENDDO 

ALPHA = 2.0D+0 


AA = 1.1D+0 


200 ALPHA = ALPHA/2. 0D+0 
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C — Calculate the parameters using the modified step-size 
DO K = 1,NP 

BS(K) = BSS(K) + ALPHA*(B(K) - BSS(K)) 

ENDDO 
CHANGE = 0 
G = 0.0D+0 

C — Calculate the slope, G 
DOK= 1, NP 
DELTAB = BS(K) - BSS(K) 

G = G + DELTAB *SUMG(K) 

ENDDO 

C — By the definition of G, it should always be positive. 

IF (G.LT.0.0D+0) THEN 

WRITE(7,*)’G IS NEGATIVE, TERMINATE CALCULATIONS’ 
GOTO 1000 
ENDIF 

SSYP = 0.0D+0 

C - Calculate the new sum of squares based on the modified parameters 
DO 13 = 1, N 
1 = 13 

CALL MODEL 

RESID = Y(I) - ETA 

SSYP = SSYP + RESIDUES ID/SIG2(I) 

ENDDO 

IF(SSYP.GT.SSY) THEN 

IF (ALPHAJLE.O.OlDfO) THEN 
WRITE(7,250)ALPHA,SSYP,SSY 

250 FORMAT(3X, ’ALPHA IS TOO SMALL, ALPHA =’,F12.6,2X, 

+ ’SSYP = E15.6, 2X, ’SSY =’, E15.6) 

GOTO 1000 
ELSE 
GOTO 200 
ENDIF 
ENDIF 

C -- Calculate SUMCH, used in the following inequality to determine H 
SUMCH = SSY - ALPHA*G*(2.0D+0-(l .0D+0/AA)) 
IF(SSYP.GT.SUMCH) THEN 

H = (ALPHA* ALPHA*G)/(S S YP - SSY + (2.0D+0*ALPHA*G)) 

ELSE 

H = ALPHA* AA 
ENDIF 

C — Calculate the final parameter estimates using H. 

DOK=l, NP 

B(K) = BSS(K) + H*(B(K) - BSS(K)) 

ENDDO 

C — Calculate RATIO; if it is less that CRTIER (0.0001), then the change 
C in the estimated parameters is insignificant and the iterative 
C process is terminated. CHANGE is used to determine when all 
C parameters stop varying. 
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DO J = 1, NP 

RATIO = (B(J) - BSS(J))/( BSS(J)+EPSDEN) 

RATIO = ABS (RATIO) 

EF(RATIO.LE.CRITER)THEN 
CHANGE = CHANGE + 1 
ENDIF 
ENDDO 

C — Print out the calculated values for H, G, SSY, and SSYP 
WRITE(7,120) 

WRITE(*,120) 

120 FORMAT(5X,’MAX’,10X,’H’,13X,’G’,12X,’SSY’,11X,’SSYP’) 

WRITE(*,125)MAX,H,G,SSY,SSYP 
WRITE(7,125)MAX,H,G,SSY,SSYP 
125 FORMAT(I8, 1F1 3 .6,4E14.6) 

C — Print out the final parameter estimates 

WRITE(7,*)’THE FINAL PARAMETER ESTIMATES FOR THIS ITERATION ARE’ 
WRITE(*,’(10X,”B(”,I1,”) =”316.6)’) (1300,1=1, NP) 

WRITE(7,’(10X,”B(”,I1,”) =”316.6)’) (I,B(I),I=1,NP) 

C — End the Box-Kanemasu Modification 

WRITE(7,’(/,5X,”P(1,KP)”,9X,”P(2,KP)”,9X,”P(3,KP)”,9X, 

+ ”P(7,KP)”,9X,”P(5,KP)”)’) 

C — Print out the P matrix 
WRITE(7,129) 

129 FORMAT(5X, ’THE P MATRIX IS’) 

DO IP = 1, NP 

WRITE(7,130) (P(IP,KP),KP= 1 ,NP) 

ENDDO 

130 FORMAT(5D15.7) 

WRITE(7,135) 

135 FORMAT(5X,”IHE CORRELATION MATRIX IS’) 

DO IR=1,NP 
DO IR2 = 1, IR 
AR = P(IR,IR) * P(IR2,IR2) 

R(IR,IR2) = P(IR,IR2)/SQRT(AR) 

ENDDO 
ENDDO 
DO IR = 1, NP 

WRITE(7, ’ (5E15 .7)’ )(R(IR,IID,ni= 1, IR) 

ENDDO 

WRITE(7,*) 

WRITE(7,*)’The diagonal terms of the correlation matrix are 
+ all unity and the off-diagonal terms must be in the interval 
+ [-1,1]. Whenever all the off-diagonal terms exceed 0.9 

+ in magnitude, the estimates are highly correlated and 
+ tend to be inaccurate’ 

DO K = 1, NP 
XTY(K) = 0.0D+0 
BS(K) = B(K) 

SUMG(K) = 0.0D+0 
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DO J = 1, NP 
XTX(J,K) = 0.0D+0 
PS(J,K)=O.0D+O 

PS(JJ) = BSV(J)*BSV(J)*1.0D+7 
ENDDO 
ENDDO 
WRITE(7,400) 

400 FORM AT(7X, ’MAX’ ,8X, ’ NP’ , 5X, TP’ ) 

WRITE(7,’(7I10,4F10.4)’)MAX, NP 
IF(NP.GT.CHANGE)T1IEN 
M = MAX IT 


IF(MAX.LE.M)GO TO 100 
ENDIF 

IF(IPRINT.LE.O) THEN 
IPRINT = IPRINT + 1 
ENDIF 

1000 CONTINUE 
CLOSE(7) 

CLOSE(8) 

STOP 

END 


CCCCCCCCC************************************************CCCCCCCCC 


SUBROUTINE MODEL 

C THIS SUBROUTINE IS TO CALCULATE T, THE TRUE TEMPERATURE 
IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION T(1500,5),Y(1500),SIG2(1500),B(5),Z(5), 
+A(5),BS(5),VINV(5,5)3XTRA(20) 

DIMENSION P(5,5),PS(5,5),SUM(5) 

COMMON T,N,Z,BS,I,ETA,PS,P,B,A,Y,SIG2,MODL,VINV,NP,EXTRA 

COMMON/MOD/AA,TL,SUM 

PI=4.0D+0*DATAN(1 .0D+0) 

QO = EXTRA(l) 

TO = EXTRA(2) 

FIMEH = EXTRA (3) 

AL = EXTRA(4) 

AL2 = AL*AL 


X = 0.0D+0 


TIME = T(I,1) 

THCON = BS(1) 

RHOCP = BS (2) *1.176 
TOL = 1J7-8 
XL = X/AL 

C DIMT = (ALPHA*t/L A 2) 

DIMT = (THCON*TIME)/(RHOCP*AL2) 

DIMTII = (THCON*(TIME-TIMEH))/(RHOCP*AL2) 
C CONST = (QL/K) 

CONST = (QO*AL)/THCON 
SUMT = 0.0D+0 
DO 20, M = 1, 1000 
BETAM = (M - 0.5D+0)*PI 
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BETA2M = B ET AM * B ET AM 
FF1 = DEXP(-BETA2M*DIMT) 

IF(TIMEJLE.TIMEH) THEN 
T1 = FF1 
ELSE 

FF2 = DEXP(-BETA2M*DIMTH) 

T1 = FF1 - FF2 
ENDIF 

TINCR = Tl*DCOS(BETAM*XL)*(l/BETA2M) 

IF(M.NE.l) THEN 

IF(ABS(TINCR/SUMT).LT.TOL) THEN 
GOTO 15 
ENDIF 
ENDIF 

SUMT = SUMT + TINCR 
20 CONTINUE 
15 IF(TIME .LE.TIMEH) THEN 

ETA = TO + CONST*(1.0D0 - XL - 2.0D0*SUMT) 

ELSE 

ETA = TO - 2.0D0*CONST*SUMT 
ENDIF 
GOTO 1000 
1000 CONTINUE 
RETURN 
END 

SUBROUTINE SENS 

C THIS SUBROUTINE IS FOR CALCULATING THE SENSITIVITY COEFFICIENTS 
IMPLICIT REAL *8 (A-H,0-Z) 

DIMENSION T(1500,5),Y(1500),SIG2(1500),B(5), 
+Z(5),A(5),BS(5),VINV(5,5),EXTRA(20) 

DIMENSION P(5,5),PS(5,5),SUM(5) 

COMMON T,N,Z,BS,LETA,PS,P,B,A,Y,SIG2,MODL,VINV,NP,EXTRA 

COMMON/MOD/AA,TL,SUM 

PI=4.0D+0*DATAN(1 .0D+0) 

TZ=0.0 

QO = EXTRA(l) 

TO = EXTRA(2) 
riMEII = EXTRA (3) 

AL = EXTRA (4) 

AL2 = AL*AL 
X = 0.0D+0 
TIME = T(I,1) 

THCON = BS(1) 

RHOCP = BS(2)*1.D6 
TOL = 10-8 
XL = X/AL 

C DIMT = (ALPHA* t/L A 2) 

DIMT = (THCON*TIME)/(RHOCP*AL2) 

DIMTH = (TIICON*(TIME-TIMEH))/(RiIOCP*AL2) 
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C CONST = (QL/K) 

CONST = (QO*AL)/THCON 
SIJMK = 0.0D+0 
SUMC = 0.0D+0 
DO 20, M = 1, 1000 
BET AM = (M - 0.5D+0)*PI 
BETA2M = BETAM*BETAM 
FF1 = DEXP(-BETA2M *DIMT) 

IF(TIME.LE.TIMEH) THEN 
XI = FF1 *((1/BETA2M) + DIMT) 

X2 = DIMT*FF1 
ELSE 

FF2 = DEXP(-BETA2M*DIMTH) 

XI = ((1/BET A2M)+DIMT)*FF1 - ((1/BET A2M)+DIMTH)*FF2 
X2 = DIMT*FF1 - DIMTH*FF2 
ENDIF 

X1INCR = XI *DCOS(BETAM*XL) 

X2INCR = X2*DCOS(BETAM*XL) 

1F(M.NE.1) THEN 

IF(ABS(XlINCR/SUMK).LT.TOL.AND.ABS(X2INCR/SUMC).LT.TOL)THEN 
GOTO 15 
ENDIF 
ENDIF 

SUMK = SUMK + X1INCR 
SUMC = SUMC + X2INCR 
20 CONTINUE 
15 EF(TIME.LE.TIMEH) THEN 

Z(l) = -(CONST/THCON)*(1.0D0 - XL - 2.0DO*SUMK) 

Z(2) = -2.0DO*(CONST/RHOCP)*SUMC*1.D6 
ELSE 

Z(l) = 2.0DO*(CONST/THCON)*SUMK 
Z(2) = -2.0DO*(CONST/RHOCP)*SUMC* 1 .D6 
ENDIF 
GOTO 1000 
1000 CONTINUE 
RETURN 
END 
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Appendix E 


Input file for MODBOX.FOR 

This file represents a sample input file to be used in MODBOX.FOR for the 
estimation of the thermal properties. The first row of numbers represents the number of 
data points, the number of parameters to be estimated, the number of independent 
variables, the maximum number of iterations to be performed, the model number, and the 
usual printouts, respectively. The second row represents the initial guesses for k x _ eff and 
C efft which are to be estimated. The first column is the index, the second column is the 
values of the temperatures, the third is the standard deviation of the measurement errors, 
and the fourth column is the independent variable, time. At the end of the columns, the 
first row represents the number of extra parameters to be used in the program. These 
parameters are then given in the next row, and represent the magnitude of the heat flux, 
the initial temperature (and in this analysis, the constant temperature at the boundary), the 
heating time, and the composite thickness, L x , respectively. 

1045 2 1 2 1 

0.5 1.5 

1 .201515E+02 .100000E+01 .OOOOOOE+OO 

2 .202875E+02 .100000E+01 .500000E+00 

3 .204508E+02 .100000E+01 .100000E+01 

4 .20547 1E+02 .100000E+01 .150000E+01 
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5 .207033E+02 .100000E+01 .200000E+01 

6 .207728E+02 .100000E+01 .250000E+01 

7 .209213E+02 .100000E+01 .300000E+01 

8 .209658E+02 .100000E+01 .350000E+01 

9 .210080E+02 .100000E+01 .400000E+01 

10 .211043E+02 .100000E+01 .450000E+01 


1035 .201841E+Q2 .100000E+01 .517086E+03 

1036 .202061E+02 .100000E+01 .517586E+03 

1037 .2021 14E+02 .100000E+01 .518086E+03 

1038 .202536E+02 .100000E+01 .518586E+03 

1039 .201520E+02 .100000E+01 .519086E+03 

1040 .20283 3E+02 .100000E+01 .519586E+03 

1041 .20241 1E+02 .100000E+01 .520086E+03 

1042 .202583E+02 .100000E+01 .520586E+03 

1043 .201989E+02 .100000E+01 .521086E+03 

1044 .202584E+02 .100000E+01 .521586E+03 

1045 .201692E+02 .100000E+01 .522086E+03 

4 

351.05 20.1515 180.036 0.0067818 
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Appendix F 


Finite Element (EAL) program 


This is the finite element program (EAL) used to estimate the thermal properties. 


$ EAL THERMAL ANALYSIS RESEARCH PROJECT, EXP. 3 

$ 

$ Debbie A. Moncman 

$ 

$******* ******************************************************* ************ 

$ 

$ This problem solves for the temperature distribution in a 2-D plate 
$ with dimensions LXm x LYm. It then uses the Modified Box-Kanemasu 

$ method to sequentially estimate the thermal properties of the material. 

$ The properties of interest are the effective thermal conductivity and the 
$ volumetric heat capacity. The left and right surfaces of the flat plate are 
$ insulated, the bottom surface is maintained at a const, temp, and at 
$ the top surface is a constant heat flux. The assumptions 
$ used are: Transient, one-dimensional conduction, constant properties, 

$ and no internal heat generation. 

$ 

£************************************************************************** 


*XQT U1 


*CM=200000 

$ 

<£***+** ****************************** ************************ ************** 


Subroutine VARB - defines variables used in the program 
NOTE: Variable names can only be four 


$ letters long! 

Q*********** *********************************))( lit******** ************** ****** 

$ 

*(29 VARB DEH) VARB 


$ 

$ Set RACM = 0 to use Fortran logic in all subroutines 
*RACM = 0 
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$ 

$ Define geometry for 2-D plate 

$ 

!LX=0.0479425 SLength of plate (m) 

!LY=0.0067818 $Height of plate (m) 

$ 

$Define number of elements and nodes in each direction 

$ 

!NX=5 $Number of elements in X direction 

!NX1=NX+1 $Number of nodes in X direction 
!NX2=NX1+1 $JJUMP for meshing; start of second row 
!NY=5 $Number of elements in Y direction 

!NY1=NY+1 $Number of nodes in Y direction 
!NT=NX1*NY1 $Total number of nodes in mesh 

!NTOT=NT 

!N1=1 $Beginning node for mat’l 1 (only 1 mat’l in this analysis) 

!RN1=1. 

!RINC=1. SNodc number increment 

!CRIT = l.E-6 $Criteria used to terminate estimation process 

!IT1 = 0 $Value used in determining if ests. are still changing signific. 

!IT2 = 0 $Value used in determining if ests. are still changing signific. 

$ 

$ Define initial temperature, initial and final time, time step for 
$ transient solution, total heating time, and heat flux value. 

$ 

!TEMI=20.1515 $Initial temperature (oC) 

! TIM 1=0.0 $Initial (starting) time (sec) 

!TIMF=522.086 $Final (stopping) time (sec) 

!DELT=0.5 $Time step for transient solution 

!TTMH=180.036 $Time that heat flux is applied. 

!DFLX = l.E-8 $Small value added to timeh to define heat flux value 

!THDL = TIMH + DFLX $This is needed due to the discontinuity at timeh 
!FLUX = 350.05 $Heat flux value 

$ 

!COU=TIMF-TIMI 

!COU=COU7DELT SNumber of time steps used in taking temp, measms. 

!NTS=1WX(COU+0.0001) SNumber of time steps must be an integer 
!NTS=NTS-fl STotal number of times from mil to TEMF 

$ 

$ Enter initial parameter estimates 

$ 

!A1I=2.0 Slnitial estimate for eff. thermal conductivity 

!A2I=3.0 Slnitial estimate for eff. volumetric heat capacity 

!LOOP=0 SUsed in sequential process 
♦RETURN 
*VARB 
$ 


^ + + + * + + + + + + + + + + + + + + * + + + + * + + + + + + + + + + + + + + 


S 


$ Subroutine NODE - defines the nodal positions 
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$ 

£**************************************************************************** 

$ 

*(29 NODE GENE)NODE 
*XQT TAB 

START "NT" $ Define the total number of nodes 
UPDATE=1 

JLOC $ Give the location of the nodes (set up the mesh) 

$ 

$ In the next statement, FORMAT=l is used for rectangular coordinates; 

$ N1 is the number to start the node locations at (in this case, 1); 

$ 0,0,0 are the coordinates of Nl; LX,0,0 are the coordinates 
$ for the bottom right comer of the mesh; NX1 is the number of nodes 
$ in the x direction; 1 is the increment in the node number in the x 
$ direction; and NY1 is the total number of nodes in the y direction. 

$ For the next line, NX1 is jjump used in the y direction; 0, LY, 0 are 
$ the coordinates of the upper left node; and LX, LY, 0 are the 
$ coordinates of the upper right node. 

$ 

FORMAT = 1: "Nl", 0., 0., 0.,"LX", 0., 0.,"NX1",1,"NY1" 

"NX1”, 0., "LY", 0.,"LX", "LY", 0. 

♦RETURN 

♦NODE 


$ 

<jj * * * * * * * * * * + * * * * * * * * * * 9|C * * ** * * * * * * * * * * * ** * * * * * * * * ]>( * * * ** « * * * * * 4c * % * * * * % * * * * * * * * 

S 

$ Subroutine ELEM - defines the element connectivity 
$ 

$ * * ******* ********** ******************************************************** 
$ 


*(29 ELEM DEFQELMT 
♦XQT ELD 
$ 

$ Define K41 elements 
$ K41 signifies a conductive, 4 node element. 

$ 

RESET NUTED=1 
K41 

GROUP = 1 $ Group 1 

NMAT=1 $Material 1 

!J1=N1 $J1 is the number of the bottom LF node in an element 

!J2=N1+1 $J2 is the number of the bottom RT node in an element 


!J3=NX2+1 $J3 is the number of the top RT node in an element 

!J4=NX2 $J4 is the number of the top LF node in an element 

"Jl" "J2” "J3" "J4", 1, "NX" "NY" 


$ 

$ 

$ 

$ 

$ 


The above line sets up the nodal positions of each element. 1 is the 
node increment and NX and NY signify the number of elements in the x and 
y directions, respectively. 
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$ Define K21 elements 
$ K21 is used to represent the heat flux 

$ Note: since we’re using this element to model the heat flux, the 
$ thermal conductivities and specific heat must be zero (reason for Mat’l 2) 

K21 

GROUP = 1 
NMAT = 2 
!J1=NY*NX1+1 
!J2=NY*NXl+2 
"J1 M ,, J2 H , 1, "NX", 1 
$ 

♦RETURN 

*ELMT 

$ 

$ 

$ Subroutine TABL - Defines the thickness of the elements 
$ 

^♦♦ft*********************************************************************** 

$ 

*(29 TABL GENE)TABL 
*XQT AUS 

TABLE(NI=1,NJ=1): K THIC: J=l: 1. $The thickness of K41 elements 
TABLE(NI=1,NJ=1): K AREA: J=l: 1. $The area of K21 elements 
$ 

♦RETURN 

♦TABL 

$ 

$ 

$ Subroutine UPDA - Updates thermal property values 
$ 

$* ******* ***** *************************** *** ****** Sic********** ********* ****** 
$ 

*(29 TABL UPDA) UPDA 
♦XQT AUS 
$ 

$ The following table gives the thermal conductivity in the x and y direc. 

$ NI=9 indicates that nine variables can be entered to determine k (T, rho, 

$ c, kxx, kyy, kzz, kxy, kzy, kzx) however, NJ=1 indicates that k is 
$ temperature independent. 1 = 45: correspond to kxx and kyy inputs (i.e. 

$ the thermal conductivities in the x and y directions). NOTE: all non-zero 
$ conductivities must be specified; there are no default values. To define 
$ isotropic properties, identical values for kxx, kyy, and kzz must be entered. 

$ Note that a value of 1.0E+6 was given for the density. This is used as a 
$ scaling factor. Therefore, the estimate for the volumetric specific heat 
$ must be multiplied by (1.0E+6). 

$ 

TABLE(NI=9,NJ=1): COND PROP 1: 1 = 2,3,4,5,6 
OPERATION=XSlJM 
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J=l: 1.E+6, ,, A2","A^7'A^ , ,"A^ , 

$ 

$ The next table sets up the 0 value properties for the K21 element 
$ 

TABLE(NI=9,NJ==1): COND PROP 2: 1 = 3 4 5 6: J=l: 0.,0.,0.,0. 

$ 

$ The next table defines the constant heat flux applied to the top surface. 

$ Here, J is the number of elements. To specify a line heat flux (W/m) along 
$ prescribed line flux divided by the element’s cross-sectional area. 

$ A constant heat flux is applied for a timeh "TIMH", and is then removed. 

$ (It is a step function) 

$ 

TABLE(NI= 1 ,NJ="NX"): SOUR K21 1: BLOCK 1: J=l, "NX": "FLUX" 

BLOCK 2: J=l, "NX": "FLUX" 

BLOCK 3: J=l, "NX”: 0.0 
BLOCK 4: J=l, "NX": 0.0 
TABLE(NI-1,NJ=4): SOUR TIME: J=l: "TIMI" 

J=2: "TIMH" 

J=3: "THDL" 

J=4: "TIMF" 

$ The next table defines the nodes that have a prescribed temperature 
$ (The bottom surface in this analysis). DDATA is a counter; Note: it must 
$ be a REAL value (not an integer). Here, J is the number of nodes, NOT 
$ the node number! 


TABLE(NI=1,NJ="NX1 "): TEMP NODE: DDATA="RINC":J="N1","NX1": "RN1" 
$ 

$ The next table defines the prescribed temperature at each of the nodes 
$ specified in the previous table. In this analysis, all specified nodes 
$ are at the same temperature "TEMI" since exps. were conducted at room temp 
$ 


TABLE(NI=1,NJ="NX1"): APPL TEMP: BLOCK 1: J=1,"NX1": "TEMI" 

$ 

♦RETURN 

♦UPDA 

$ 

^********* *********************************************************** ******* 
$ 

$ Subroutine TDAT: Builds data files for experimental temperatures, initial 
$ guesses, and sensitivity coefficients. 

$ Data format: 

$ TS AUS n3 1 

$ Specify both n3 and n4 to identify the data tables 

$ n3 = 1: measured temperatures 

$ n3 = 2: initial temperatures 

$ n3 = 3: derivative #n (the n-2th parameter) 

$ n4 - x position where temperatures are measured 

$ (only measured at one location (n4=l) for this case). 

^** ***************************************************************** ******** 
$ 
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*(29 BILD TS) XXXX 
1NC02 = NCOU 
$ 

$ Bring surface temperature data from TRAN TEMP multi-block data set using 
$ XSUM and TRANSFER. 

S 


£ * 9|c * ** * * * * ** * * * * * * * * * * * * * * * * **** ******** * * * ***** * * * * *** *** * ** >|e ** * * ** * * ** s|c * * * 
$ 


!NTN = NTOT-1 $Total number of nodes - 1. 

!NCP = NCOU+1 $=NTS, the total number of times from TIMI to TIME 
!DBS = 0 


*IF("N4" EQ 2):!DBS=NCOU 
*XQT AUS 

DEFINE A = 1 TRAN TEMP 1 1 2 "NCP" 

*JFCNT GT 1):*JUMP 1215 

*IFrN4 H EQ 1):TABLE(NI=1,NJ="NCP , '): 4 TS AUS "N3" 1 
*fF("N4" EQ 2):TABLE,U (NI=1,NJ=’ , NCP ,, ):4 TS AUS "N3" 1 
TRANSFER(SOURCE=A r DBASE="DBS , ',SBASE="NTN",ILlM=l,OPERATION=XSUM) 
♦LABEL 1215 

*IF("N4" EQ 1):TABLE(NI=1,NJ="NCP"): 2 TS AUS "N3" 1 
♦EFO'W EQ 2): TABLE, U (NI=1,NJ="NCP"): 2 TS AUS "N3" 1 
TRANSFER(SOURCE=AJDBASE= M DBS ,, ,SBASE="NTNMLlM=l,OPERATION=XSUM) 
♦RETURN 


♦XXXX 


* * * * * * * * * * ♦ * ♦ * * * s|t * £ * * ♦ ♦ ♦ * * * * ♦ * * ♦ * * * * * * * * * * * * * * * sfc * iff lie $ * * * * afe a|t sft * * jfcafc * * * * * * * £ >fc * >|c * 

$ 

$ Subroutine INVH - Minimization Procedure 

$ 

^********* ************************************************ ****************** 

$ 

*(29 INV HEAT) INVH 
$ 

!NCOU=NTS-l $Total number of time steps 


♦XQT AUS 
$ 

TAB LE(NI= 1 ,NJ= 1045): 4 TS AUS 1 1: 1=1 


J=l:.201515E+02 
J=2: .202875E+02 


J=3:.204508E+02 


J= 1043 : .201989E+02 
J=1044:.202584E+02 
J=1045: .201692E+02 
!NCOU=NTS-l 

!EPS=1.0E-6 $Convergence criteria used in (b-b0)/(b0-EPS) 

!NEPS=1 $Used to determine is No. of iterations exceeds NEMX. 
!NEMX=25 $Maximum number of iterations 
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*XQT AUS 

!TINI=DS 1,1,1 (4,TS AUS, 1,1) $Defines the initial exp. temperature 
TABLE (NI=4,NJ="NEMX"): 4 CONV HIST 1 1 $Table that stores sequential est. 
STABLE (NI= 1 ,NJ="NCOU"): 4 RES HIST 1 1 
!A1=A1I SSet A1 and A2 back to the initial estimates. 

!A2=A2I 

!AS1 = All $AS1 and AS2 are previous iteration, final estimate holders. 

!AS2 = A2I 

♦LABEL 4000 SBegins the loop process 
!A10=A1 
!A20=A2 
$ 

SDerivative calculations (Used to calculate the sensitivity coefficients 
1TINI = DS 1,1,1 (4,TS AUS, 1,1) 

!N4=1 

!N3=2 

!NTAB=0 

♦DCALL(29 TRAN ANAL) 

$The above call stmt calculates temps, at the initial parameter estimates for 
Sthe first iteration and at the final estimates of the previous iteration for 
Sthe 2nd, 3rd, ... NEMX iterations. 

!A1=1.001*A10 SEstimate at Al+dAl 

!DA1=0.001*A10 SStep used to numerically differentiate 

!N3=3 

♦DCALL(29 TRAN ANAL) SCalculates temps, at (Al+dAl) 

!A1=A10 SSet A1 back to initial estimate 

$ 

!A2=1.001*A20 SEstimate at A2+dA2 

!D A2=0.001 * A20 

!N3=4 


*DCALL(29 TRAN ANAL) SCalculates temps, at (A2+dA2) 
!A2=A20 SSet A2 back to initial estimate 

$ 

$ **** INVERSE HEAT TRANSFER BEGINS HERE *♦♦♦ 

$ 

$ The parameters are initially estimated using the Gauss Method. These 
$ estimated values are then modified using the Box-Kanemasu Method. 

$ 


♦XQT AUS 

INLIB = 2 Sldentifies the source library 

OUTLIB as 2 Sldentifies the destination library for output datasets 


DEFINE TM = 4 TS AUS 1 1 
DEFINE TO = 2 TS AUS 2 1 
DEFINE TA1 = 2 TS AUS 3 1 
DEFINE TA2 = 2 TS AUS 4 1 
$ 


SExperimental Temperatures (Y) 

SCalc. temps, at initial parameter est. (ETA) 
STemps at A1 + DAI 
STemps at A2 + D2A 


!A1NV = 1.0/DAI 
1A1N2 = -1.0/DAI 
!A2NV = 1.0/DA2 
IA2N2 = -1.0/D A2 


SThe following 4 statements are used in finding 
Sthe sensitivity coefficients. (The derivative of 
Sthe temperature with respect to the parameter). 
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$ 

D1 = SUM ("A1NV" TA1, "A1N2” TO) $Delta Ti (n@Al+DAl-TI@Al) 

$ This statement sums the derivatives for the thermal conductivity 
$Meaning: A1NV * TA1 + A1N2 * TO 

D2 = SUM ("A2NV” TA2, "A2N2" TO) $Delta Ti (n@A2+DA2-TI@A2) 

$This statement sums the derivatives for the volumetric heat capacity 

$ 

N1 = SUM (TM,-1. TO) $Gives the matrix (Y - ETA); the Residuals 

$ Build up the X matrix using vectors containing 

$ the derivatives 

$ 

SENS MATRIX = UNI0N(D1JD2) $Joins D1 and D2 into a new dataset 
$ D1 and D2 must have the same block length. 

$ 

DEFINE S=SENS MATRIX 1 1 SDefines the matrix X, i.e., the Sens. Coeffs. 
ERR = XTY(S,N1) $Calculates XT (Y - ETA) 

STS = XTY(S,S) $Calculates XIX 

STSI = RINV(STS) $ Calculates the INVERSE of (XTX) 

DA = RPROD(STSI,ERR) $Calculates INV(XTX)*(XT)*(Y - ETA) 

$ 

NTN = XTY(N1,N1) $Calcs the Sum of Squares, (Y-ETA)T(Y-ETA) 

TIT = XTY(TM,TM) $Calculates YTY 

!DA1 = DS 1,1,1 (2, DA AUS, 1,1) 

$DA1 is the perturbation for the new estimate (thermal conductivity) 

!DA2 = DS 2,1,1 (2, DA AUS, 1,1) 

$DA2 is the perturbation for the new estimate (volumetric heat capacity) 

!SYS = DS 1,1,1 (2, NTN AUS, 1, 1) $The sum of squares value 

$ 

$The following (A1 & A2) are the estimates obtained with only the Gauss Method 

$ 

!A1 = DA1+A1 $New parameter estimate for the thermal conductivity 

!A2 = DA2+A2 $New parameter estimate for the volumetric heat capacity 

$ 

$ *** END BASIC LOOP-BEGIN BOX-KANEMASU MODIFICATION *♦* 

$ 

$ This section of the program takes the estimated parameter values found 
$ using the Gauss Method and modifies them using the Box-Kanemasu method. 

$ This method may allow for convergence of the parameters when the Gauss 
$ method does not. It uses the direction provided by the Gauss method but 
$ modifies the step size by introducing a scalar interpolation factor (H). 

$ The final parameter values are calculated using the Box-Kanemasu method. 

$ For a detailed explanation of this method, see ’Parameter Estimation’ by 
$ J. Beck and K. Arnold (p. 362-367). 

!AG1 = A1 $Fixes the Gauss estimates 

!AG2 = A2 

!ASS1 = AS1 $Fixes the initial estimate for that iteration 

!ASS2 = AS2 

!ALPH = 2.0 $Used in finding the parameter estimates 

!AA = 1.1 $Used to calculate H 

♦LABEL 620 
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*IF("ALPH" LT 0.01):*JUMP 4001 $Alpha is too small, estms. aren’t converging. 
!ALPH = ALPII/2.0 $Alpha starts out as 1.0. 

!DIF1 = A1 - ASS1 $Diff btw. Gauss & final est. of previous iteration. 

!DIF2 = A2 - ASS2 

!AS1 = ASS1 + ALPH*DIF1 $Est using the modified step-size 

!AS2 = ASS2 + ALPH*DIF2 

!ALPHA 

DATR = RTRAN(DA) $Transpose (XTX) A (-1)XT(Y - ETA) 

G = RPROD(DATR,ERR) $Used in calc. H, it’s the slope of the Sum of Squares 
$ vs. H. By defn., it should always be a positive scalar 

1GVAL = DS 1,1,1 (2, G AUS, 1, 1) $Gives the scalar value found for G 
!A1 = AS1 
!A2 = AS2 
!N3 = 5 
!N4 = 1 

*DCALL(29 TRAN ANAL) 

$The above call stmt calculates the temp, at the estimates obtained using ALPH 
*XQT AUS 
INLIB = 2 
OUTLIB = 2 
S EXIT 

DEFINE TOG = 2 TS AUS 5 1 $Temperatures at the Gauss est. + Step Size(alph) 

DEFINE TM = 4 TS AUS 1 1 SExperimental temperatures 

NSS = SUM(TM,-1. TOG) $New (Y-ETA) using TOG temperatures. 

SYP = XTY(NSS,NSS) $New sum of squares 

!SSYP = DS 1,1,1 (2, SYP AUS, 1, 1) $Gives the sum of squares value 
$ 

*IF("SSYP" GT "SYS"):*JUMP 620 

$ 

$ The above statement is a check to see if the sum of squares is decreasing 
$ If it’s not, alpha is decreased by 1/2. This process continues until the 
$ above if statement is no longer true or until alpha is < 0.01, in which 
$ case the program is terminated. 

$ 

!CHEK = SYS - ALPH*GVAL*(2.0 - (1/AA)) $This is a check used to determine H 
!H = ALPH*AA $Mtially set the step-size, H equal to alpha* AA. 

$ 

$If SSYP > CHEK, H is given a new value; see following IF stmt. 

$ 

*IF("SSYP" GT "CHEK"):!H = (ALPH*ALPH*GVAL)/(SSYP-SYS+(2.0*ALPH*GVAL)) 
$ 

$Calculate the modified parameter values using the obtained step-size (H). 

!A1 = ASS1 + H*(AG1 - ASS1) $Parameter estimates obtained using B-K method. 

!A2 = ASS2 + H*(AG2 - ASS2) 

$ 

SCalculate the following ratios, if RATI and RAT2 are < CRTT (0.0001), then 
$the change in the estimated parameters is insignificant and the iterative 
Sprocess is terminated. 

$ 

IRAT1 = (A1 - ASS1)/(ASS1 + EPS) 
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!RAT1 = ABS(RATl) 

!RAT2 = (A2 - ASS2)/(ASS2 + EPS) 

!RAT2 = ABS(RAT2) 

$ 

!LOOP=LOOP+l $Next iteration 
♦XQT AUS 
$ 

$ Updates the table of the sequential estimates 
TABLE,U(TYPE=-2): 4 CONV IIIST 1 1 
J=”LOOP": "Ar , ,"A2","AGr , , ,, AG2" 

$ 

$Set this iterations final estimates equal to the initial estimates for 
$the next iteration. 

!AS1 = A1 
!AS2 = A2 
$ 

♦IFfRATl" LE ”CRIT"):!m = 1 
*IF("RAT2" LE "CRrF):!IT2 = 1 

!ITER = IT1 + IT2 $ Determines if the change in both ests. is insignf 

*IF("ITER" EQ 2):* JUMP 4001 $If ests. no longer change, stop iterating. 

!NEPS=NEPS+1 $Goes to next iteration 

♦IFC'NEPS" GE "NEMX"):*JUMP 4001 

$If the parameters don’t converge before the max. No. of iters., end process 

♦XQT DCU 

PRINT 2 TS AUS 2 1 

PRINT 2 N1 AUS 1 1 $Prints out the residuals for each iteration. 

PACK 1 
ERASE 2 

♦JUMP 4000 $Est. haven’t converged yet, go to next iteration 

♦LABEL 4001 $To end iteration process 

♦XQT DCU 
PRINT 4 TS AUS 1 
PRINT 4 CONV HIST 1 1 
$ PRINT 1 TRAN TEMP 1 1 
$ PRINT 1 TRAN TIME 1 1 

$The above libraries are only printed for the final iteration. (4 TS AUS 1 
$is the for each iteration; experimental temperatures). 

♦RETURN 

* INVH 
$ 

^fr************************************************************************** 

$ 

$ Subroutine TRAN - Solves direct problem using TRTB processor 
$ 

$ 

♦(29 TRAN ANAL)TRAN 

* DC ALL (29 NODE GENE) S Generate the nodes used in the mesh 
*DCALL(29 TABL GENE) $ Generate tables needed in analysis 

* DC ALL (29 TABL UPDA) $Update the thermal properties (estimates) 


195 



*DCALL(29 ELEM DEFI) SDefines the elements (Cond.,Conv., Heat Source, etc 

$ 

*XQT TGEO $Element geometry processor; it computes local coordinates 

$ and performs element geometry checks. The user MUST 

$ execute TGEO after each execution of ELD. 

♦XQT TRTB $Transient analysis processor - Implicit with C.N. code 

RESET PTV=0.00001 Tl="TTMr T2= ,, 'nMF" DT= "DELT" PRINT=0 MXNDT=1 00000 
TEMP="TEMT 
TSAVE="DELT" 

$ 

♦XQT AUS 

!NCOU=NTS-l 

!NBLO=l 

♦DCALL (29 BILD TS) 

$ 

♦XQT DCU $Processor that performs an array of database utility 

$ functions (see Manual, Section 12-1) 

DISABLE 1 EKS B 
♦RETURN 
♦TRAN 
$ 

^******* ************************************************************ ******** 

$ 

$ Main program 

$ 

^* ******************************************************************** ****** 

$ 


♦DC ALL (29 VARB DEFI) 
♦DC ALL (29 INV HEAT) 
♦XQT EXIT 
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APPENDIX G 


Uncertainty Due to Experimental Measurements 


The uncertainty in the estimated thermal properties due to experimental 
measurement errors can be found from 


dR = 






dR 
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(G.l) 


where dR is the uncertainty in the thermal property being analyzed, 8X* is the uncertainty 
in the experimental variable, and the partial derivative of R with respect to X i is the 
sensitivity coefficient with respect to the measurement, X f . In the experiments conducted 
in this investigation, error could be associated with the temperature (X r ), voltage (X v ), or 
current (X 7 ) measurements. Therefore, the uncertainty in the effective thermal 
conductivity perpendicular to the fibers (k x . e ^) would be given by 
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Using this equation, a dk Xmeff of ±0.035 W/m°C was calculated. In Chapter 5, the mean 
value for k x . eff was estimated as 0.518 ± 0.028 W/m°C. This uncertainty of ±0.028 
W/m°C, associated with the 95% confidence region, is approximately 20% smaller than 
dk x _ eJT> found from the measurement errors. This result implies that the actual error 
associated with k x . eff may be larger than estimated. 
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